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On the Iterative Solution of the Matrix 
Equation AX’—J=0 


By Pentti Laasonen 


Introduction. The general eigenvalue problem 
(AA — B)x = 0 


with square matrices A and B often appears in numerical analysis, and -it is 
frequent that A and B are symmetric and A positive definite, so that the problem 
is equivalent to the following one: 


(AI — C)y = 0, 
where 

C = A+BA-4, 

y = Aix. 


It is clear that C is symmetric. However, in order to avoid the solution of two 
special eigenvalue problems, it is desirable to have a way of forming the trans- 
forming matrix A+ in terms of the known matrix A without a previous deter- 
mination of the eigenvalues and eigenvectors of A. Indeed, this is possible by an 
iterative process, based on Newton’s algorithm, for any matrix A with real 
positive eigenvalues. 
Newton’s method. Two simple algorithms can be presented for the iterative 
solution of the scalar equation 
ax? —1=0 
with a positive a; these are 
tiu1 = $x; a jax,*, 
and 
1 


Ninn = 3x5 + —] 
‘ 


Since the first algorithm does not always converge, the second one, which always 
converges for a nonvanishing initial value, will be suitably extended to a par- 
ticular class of matrices. The conditions under which such extension can be 
accomplished are described in the following theorem. 


THEOREM. Let A denote a real square matrix with real, positive eigenvalues. 
Then the algorithm 


(1) xX = BR 
XH) = 4X a 4(AX®)—1 
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110 ON THE ITERATIVE SOLUTION OF THE MATRIX EQUATION 


where k is a non-zero constant, generates a sequence of matrices which converges to 
that solution of 


(2) AX? —I=0 


which has positive eigenvalues. Moreover, the rate of convergence is quadratic. 
Preliminary remarks. Some simple properties of a class J of Jacobi matrices 
will be useful in establishing the results of this paper. Such matrices are of the form 


C * 


@o @ a@2°*** 
0 @o a&1*** Be-1 


0 O ao--- aes 


(3) 


_ [ao, @i1,@2, ***, a, |. 








is 8 0 ao | 


If two such matrices (denoted by bracket expressions as above) are multiplied, 
then 


(4) [ao, @1,@2, ***, as |° [Bo, Bi, Bo, ve fe B | 
= [aoBo, ao81 + a180, aoB2 + 0181 + a2Bo, --+, ao8s + arBs-1 + *** + Bo], 


so the class J is closed under the multiplication; moreover, multiplication of 
matrices of this class is commutative. Finally, the inverse of a matrix in J, if it 
exists, is again in J. 

The following simple result will also be useful. 


Lemma. Let A 
fa}, {8}, fe} (é =0,1,2, ---) 


denote three sequences of real numbers. Then if they are related by the relationship 
att) = gl) 4+ BO, 


and if the sequences {8B} and {¢} have limits B and 0, respectively, the sequence 
{a} converges to the limit B. 

The proof is simple, consisting first of a finite estimate for the upper bound of 
a, whereafter the final conclusion is evident. 

Finally, this section concludes with the well known Jordan matrix theorem. 
This theorem states that any square matrix A can be transformed by a suitable 
non-singular matrix U into the Jordan normal form 
™ ‘ 

Ae 


(5) UAUt=D= 








L Am 


4 


where all elements outside the diagonal submatrices A, are zero and these sub- 
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matrices have the form 


i-_ & 6 0 0 | 

=. 3s 0 0 

00 », 0 0 
A, = 

0 0 0 2 

10 0 0 0 AJ 








Convergence of the algorithm. Define the matrix Y“ by 
yo = UX®U4, 
Evidently the algorithm for Y“:s is defined by the formula 
yoo = sy + (DY), 
Now, suppose that one of the matrices Y“, say Y, has the special Jacobi form 
fry } 
rr” 


(6) yo) = vets : 








(r) 
ra” | 


~ 


where the square submatrices I’, along the diagonal have the same orders as the 
corresponding submatrices A, and, moreover, belong to the class J defined by (3). 
It is then obvious that this particular character will be maintained by all suc- 
cessive matrices Y+»), Y‘+®), .--, and, moreover, the corresponding submatrices 
Ir, are related by 


(7) 1, = 9r,© + 3,7,0)* (a = 1,2,---,m;¢ = 7,7 +1, ---). 
Now, transform this equation into the form 
(8) ATO (24+) — PM) = J, 
where the subscripts have been deleted. Then use the notations 
TO = Lao’, ary, ara, «++, ae’, 
A = [A, 1,0, ---,0], 


substitute into the equation (8) and apply the multiplication rule (4) in order to 
obtain, for the elements a; and a,‘*+”, the system of equations 


(9) hao’ (Zao — ao’) = 1, 


k k~1 k k-1 
(10) 20D aMaSt? +25 aMaft?, -—rA LaMaf, —  aMaf,. = 0 
0 r=0 0 0 


(k = 1,2, ---,5). 
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The equation (9) gives for ao‘ the formula 


1 
2rao’® 





apt) = Sao 


Hence, if the eigenvalue ao of I is positive, then the eigenvalues 
aot aot), --- converge to the positive value \~?. The equations (10) may be 
used for the successive computation of the elements a;“*”, ao“), ---, a“), 
Indeed, if one solves the equation (10) for the element a;‘+”, one obtains the 
expression 


@+1) 1 
44 ’ , , : 144 ve 
aft y=f{1- @ af’ +3 P(ag’, at”, a af, ; ay’ / eae af’ 7% 
) () 





where P is a second order polynomial of its arguments. 

An appeal to the lemma mentioned above justifies the conclusion that the 
elements a,‘ tend to a definite limit, as 7 tends to infinity. It is easy to find that 
the limit values are 


Hen agi? = (~1) oe 


ares mE 


This implies the conclusion that the sequence of the matrices Y; is convergent 
and the limit matrix has the form 


fT, 
rT, 








\ Pn J 


where the diagonal submatrices [', are Jacobi matrices of the form (3). The 
equation 
¥=4¥+4(D¥) 


is satisfied by this matrix, and therefore this is the solution of 
DT? «af. 
Similarly, the matrices X“ tend to a limit matrix X whose transform by U is Y: 
Y = UXU-. 
This matrix X is obviously the solution of 


(2) AX? —I=0. 


The assumptions made in this proof (namely, first, that one of the matrices 
Y has the particular form (6); and second, that the eigenvalues of this matrix 
are positive) are fulfilled by choosing Xo to be equal to J, since in this case also 
Yo is equal to kJ. Of course this diagonal matrix fulfills the requirements. 
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From (7) one easily derives, by using the commutative law for matrices of 
the class J, the equation 


r,t) — At = 4(7,)-71(7, — A>). 
This equation implies the further equation 


yr) — P- 


H(¥O)-(YO — Dy, 
or, finally, 
XG) — 444 


4(X@)-1(X@ — 4-4)? 


This relationship indicates that the rate of convergerce of the approximations 
XM, X®, X®, --- to the solution X is quadratic. 

The effect of round-off errors. For the quadratic rate of convergence or for 
the convergence at all of the applied Newton’s method, it is essential that the 
matrices involved have a similar structure in the sense described above; that is, 
their transforms by the same matrix A have a similar Jacobi structure (6). A 
simple counter example suffices to show that the algorithm may diverge if applied 
to matrices not related in this sense. In fact, this happens even if the initial matrix 
Xo is arbitrarily close to the correct solution. 

Such an example is provided by the matrices 


1 0 eo 
as 7. " 
. (0 ay : (6 ») 


Now, if the initial matrix is taken to be 


xo = (j a (e ¥ 0) 


then the algorithm leads to the successive approximations 


xX® = 1 (=*)*. ie 
= 2 4 (¢ = 1,2, ---). 


0 m 





This sequence diverges, if m > 3, despite the arbitrary closeness of X to the 
correct solution. In general, it can be shown that in this sense the algorithm (1), 
if carried out indefinitely, is not stable whenever the ratio of the largest to the 
smallest eigenvalue of A exceeds the value 9. 

This observation raises the question of convergence of the above method, if 
it is applied numerically with truncated values. In this case, although one starts 
with an admissible initial matrix X = kJ, the method will in general lead to a 
sequence of approximations which do not possess the required character, due to 
round-off errors. Therefore, if the method is continued indefinitely, it will probably 
not converge. For this reason it may be appropriate to compare during the process 
two successive approximations X“ and X‘*» and discontinue the algorithm as 
soon as their difference no longer decreases (by some proper norm). Due to the 
quadratic rate of convergence, in most cases this finite process does not involve 
many steps; thus, the influence of the round-off errors is likely not very great. 
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This is true, in particular, if the numerical calculations have been made with some 
additional accuracy, which can be assumed to take care of these truncation errors. 
However, it seems desirable to have a method for checking the correctness of 
an approximation and, if necessary, to improve its accuracy. This is, in fact, 
possible by the following means. 
Let X be an approximation for the exact solution X = A and form the error 


A = XAX —- I. 
If 6, 


(11) §6= A+ — X, 


denotes the correction which has to be added to X, then A can be expressed in 
terms of § and X as follows. 


A= (X +8)? = (+ X48) 7X“ fY 
= (X71 — X-9X1 + _~---)? = X* — KX — R-5X +-_—-- -. 
Thus, 
A= — X¥-§ — 6X7 +.-.--. 


The series expansions, which are convergent provided 4 is so small that the eigen- 
values of X—4 are less than unity in absolute value, are to be truncated after the 
linear terms in 6, the error caused thereby being of second order with respect to 6. 
An approximation of 6 within this accuracy may thus be obtained by solving 
the equation 


(12) X-5 + 6X99 = —A=1- XA. 


Accordingly, the succeeding paragraph is devoted to a study of the solvability 
and practical methods for solving this equation. 
The equation CV + VC = K. 


THEOREM. Let C and K be two (n X n)-matrices, all eigenvalues of C being 
different from zero and of the same sign. Then the equation 
(13) CV+ VC=K 


has one and only one solution V. 
In order to prove this theorem, transform all matrices involved by a non- 
singular matrix U, 


D = UCU, 
L = UKU-, 
W = UVU-, 


such that D = (d;;) attains any Jacobi form, that is, all elements below the 
diagonal, d; ; with i < j, vanish; the diagonal elements d;,; thus obtained are the 
eigenvalues of C and therefore of same sign, say positive. Accordingly, the equa- 
tion is equivalent to the equation 


DW+WD=L. 
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By setting corresponding elements on both sides of this equation equal to one 
another and using a self-explanatory notation, one obtains the following system 
of n? equations: 


n j-1 
4, i5+ LD dgvwes + LD werdes + weds = hs (i,j = 1,2, ---, 0). 
k=1 


kmitl 


Since the coefficient d; ; + d;,; of w;,; is positive, this equation can be considered 
as a recurrence formula which determines w,; ; in terms of the elements on the left 
lower side of the line through w;,; and parallel to the diagonal. Accordingly, all 
elements can be determined successively, beginning from the lower left hand 
corner. 

Of course, the numerical solution of the equation (13) can be achieved either 
by solving, for the elements of V, the linear system which is obtained by setting 
corresponding elements on both sides of (13) equal to one another, or by first 
performing the transformations used in the above proof. However, the first method 
is quite laborious for an arbitrary C, since the number of simultaneous equations 
involved is *?. On the other hand, the second method involves actually the 
diagonalization of the matrix C. 

A third method, based on an iteration process, is certainly preferable. The 
equation (13) is then written as follows: 


1 3 1 1 
(14) (c+tv)-c+ix+ity, 
™m™ ™m™ m 


where m is any scaling factor. If m is large enough, then a good approximation 


4 
for C +o V is provided by the matrix (c + +x) . Subtract C from this 


matrix in order to obtain an approximation for — V. This can now, if necessary, 
m 


be substituted on the right hand side of (14) and the process continued. 

The extraction of the square-root of the matrix on the right hand side of (14) 
can, of course, be achieved by Newton’s method in a way similar to the one used 
for computing the matrix A}. For the equation 


xX*-A=0, 


where A is a matrix with non-negative eigenvalues, that solution X which has no 
negative eigenvalues is obtained as the limit matrix from the algorithm 


(15) Xinn = §X; + 4AXSH, 


if the first matrix is Xp = kJ with any positive k. 

The determination of the correction 5. In applying the method outlined above 
to the determination of the correction (11) from the equation (12), observe that 
the neglected terms are of the order [(X-8)*]. Accordingly, the truncated equa- 
tion (12) should properly be replaced by 


X-§ + 6X41 = — A + Of (X8)*] 
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and therefore the corresponding equation (14) is actually 
i} 1 1 1 
Io + —38) = X* —-—A+—0[(2-9)*] + —F. 
m m m m 


Now, since the third term of the right hand side of this equation is to be 
neglected, one concludes that the removal of the last term does not cause any 
additional inaccuracy, if m is large enough. It is sufficient to choose m of the same 
order of magnitude as the largest eigenvalue of X squared or the reciprocal of 
the smallest eigenvalue of A. Compute the matrix (x _ - a) by means of 
the algorithm (15); subtract X— from this result; and finally, multiply by m in 
order to obtain an approximation for the correction (11). Of course, this pro- 
cedure may be repeated and thus set up a quadratically convergent algorithm 
which, moreover, is self-correcting. 

Summary. It has been proved that for any real (m X m)-matrix A with only 
positive eigenvalues the algorithm (1), with an initial matrix X = kJ, converges 
quadratically to the matrix A- with positive eigenvalues. 

In a numerical case this algorithm, if continued indefinitely, may be divergent 
due to round-off errors, whose influence may increase in geometrical progression. 
This makes it necessary to stop the process as soon as the difference between 
two successive results no longer decreases; of course, it is also desirable to have 
some additional accuracy in the numerical computation to take care of the 
round-off errors. 

Any approximation to A- sufficiently accurate can be improved successively, 
the rate of convergence of the procedure being quadratic. Each step, however, 
involves either the solution of a system of m? linear equations or the extraction 
of the square-root of a matrix, which may be achieved by a quadratically con- 
vergent iteration process. 


Kasarmikatu 2B11 
Helsinki, Finland 


On Gauss’ Speeding Up Device in the Theory of 
Single Step Iteration 


By Alexander M. Ostrowski 


1. In this paper we consider throughout rea] numbers, vectors and matrices. 
In order to solve the linear system 


(1) b, Quite =Jn Gee ™ A, (um = i. a -, 2) 
y=] 
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with the matrix A, A,~0 (u = 1, ---, ”) and non-vanishing determinant, by the 
single step iteration we form, starting from an arbitrary (row) vector £, a sequence 
of vectors 


(2) & = (x1, -++, Xn) (k es a 2, s+) 

obtained in the following way: For any integer k(k > 0) choose a value N; of 
the “leading index’’ from the indices 1, ---, 2; then if 

(3) pe = (rr, +++, 70) (ke = 0,1, 2, =) 


is the k-th residual vector defined by 


(4) = Fax” —y, (u =1,---,m;k =0,1,---), 
v=] 
we put 
k) 
r, + 
(5) lk +1) i xt) (u * Ni), xis) = x) Ls os. 
Ni 


In studying this iteration we can and will restrict ourselves to the case where 
all y, vanish, since we can always by a convenient change of the origin make all 
y, = 0, without changing the px. 

2. We consider in what follows only the case where the matrix of the system 
(1) is symmetric and the quadratic form 


(6) K(é) = Lo Qyx x, = EAz’ 
Byr=l 
defined for an arbitrary vector = (x1, ---, xn), is positive definite. In this case 


it is well known and immediately verified that if the vector £4, is obtained from 
the vector & by the transformation (5), we have 


(7) K (&e41) = K(&) — ()*/Aw,- 


In using (7) it was proved by Seidel, 1874 [13], that the single step iteration is 
always convergent if NV; is chosen at each step so that 


(8) (rNi2)*/Am, = maxy(rh")*/Ay 


This is Seidel’s relaxation procedure. 

This special rule goes back to F. R. Helmert, 1872 [7]. The relaxation rule 
indicated previously by Gauss [4] and Gerling [6] is different, as is the one pro- 
posed by Southwell [14], but’ the rule (8) is apparently the most advantageous 
one. Cf. the discussion in [9], p. 158-9. 

On the other hand Schmeidler [12] and Reich [11] proved in using (7) that 
the single step procedure is convergent in the cyclic case when N, runs periodically 
through all indices 1, ---, m. 

3. Gauss [4], [1], [6] and [15] proposed the following modification of the 
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above procedure in order to speed up the convergence. Put 


n 


Xy = 2% — 2, Gy =~ A0= — Lady (vy =1,---,n) 
p=l 
(9) 
n n 
ao = Ao= — Lan = DY ay, 
y=] Bel 


where 29 can be arbitrarily chosen. Then the system (1) can be written in the 
form (assuming y, = 0) 


(10) Tame (u = 0, a +++, Mm), 
oad 


where the first equation is, of course, not independent of the last equations but 
is useful for the sake of uniformity and for checking purposes. 

In particular Ao is positive since by (9) Ao is the value of the quadratic form 
(6) for x,» = 1 (v = 1, ---,m). 

4. From a solution (20, 21, -*-, Zn) of the system (10) we obtain at once by (9) 
the solution (x;, ---, 2a) of the system (1). The idea of Gauss is now to apply 
the procedure described in (4) and (5) to the system (10). If we obtain then, 
starting from a vector fo = (z0, 2:1, ---, 2,), a sequence of vectors 


ie = (zo, 2), tee, Zn), 
we consider at the same time the corresponding vectors 


fe = (25 — 29, «++, ax — zo). 


If then in the passage from ¢; to {%41 the leading index NV; ¥ 0, we have 


n n 
k k k k 
(11) x a,x” _ Y a,,(2” — 2”) _ 7") (u = 1, 7, n), 
y=0 p=] 
x k 
2 ayy,% 
(+1) _ oft) _ am _ 
zZ = zZ =— ———-— = 2N — 
; , Ay, . Ay,’ 


k+1 k 
get) = f(y x Ny). 


Since here 2**!) = 2{*), we see that the corresponding n-dimensional vectors &, 
£41 are connected exactly by the formulae (5), so that in this case there is no 
essential change compared with the original method. 

5. If, however, V;, = 0, then only zo is changed and therefore al] components 
X1, ‘**, X, are changed by the same amount. In this case we have obviously a new 
possibility and the question arises whether in this case the convergence is indeed 
speeded up. Of course, by ‘convergence’ in this case is not meant the convergence 
of the vectors {; but that of the corresponding vectors &. This question is ap- 
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parently not as yet settled, as widely contradictory opinions are to be found in 
the literature; see [14] and [3]. 

Our conclusions are stated at the start of section 8 and at the end of the 
paper. 

6. In what follows we will say that the two (m + 1)-dimensional vectors 


£ = (Zo, 21, °°, Zn) and ¢” = (z9’, ;', ---, Zn’) are equivalent, if we have z, — zo 
= 2,’ — zo’. In the class of vectors equivalent to ¢ there exists a reduced one: 
= = (0, x1, ---, xn), and the corresponding n-dimensional vector £ = (x, ---, Xn) 


is uniquely determined. 
If we define the component of the residual vector corresponding to the 
index 0 by 


n n 
(12) ro*) = > Ao,2, = 5 dy (2, = Zo) = — > 7, 
y=0 y=] y=] 


we see from (11) and (12) that the residual vector for the system (10) does not 
depend on the component 29 but only on the corresponding vector,é. It follows 
then from (9) and (6) 


n n 
X Gw2,2% = > 2, ( LD Gyre, + a,0%) 
p=0 


py y=] 


n n n n 
_ LD & DL Aux, = Lt L Gut, 


u=0 y=] y=] u=0 


= > Xy ( > Oy, + ass) 


y=l p= 


n n 
van p> Xy >> DprX yp, 


v=1 w= 


(13) p DyryZ 2, = > DyyX pXp- 


py By=l 


7. It is obvious that the algebraic identity corresponding to (7) remains true 
for the system (10), although the corresponding quadratic form is only semi- 
definite. Therefore from (13) it follows that the relation (7) is also true for Ni; = 0 
where ro is given by (12), and the quadratic form K is the positive definite 
quadratic form (6). But then it follows from (7) that in any case 


(14) lim (r$)?/Ay, = 0. 


ko 


8. We discuss first the cyclic one step iteration. In this case{we will prove 
that the procedure remains convergent, but for any m > 2 there exist matrices 
for which the modified procedure is slower and others for which the modified 
procedure is indeed faster than the original one. This agrees with results men- 
tioned by Forsythe and Motzkin [3], footnote 24. 
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9. It follows from (14) that lim r¥) = 0. Therefore from (5) and the corre- 
sponding formula for ¢ with NV; = 0, 


x) —x,© 350 (ko wo; =1,---,n), 
and by (4) 

rk) — 7.) +0 (kR- o©;p=1,---,m). 
More generally for each constant integer 7, 

rk) — 7, +0 (k— ©;p=1,---,m). 


But for any fixed u, among the m + 1 consecutive values of & there is one for 
which N; = u; therefore it follows that 


O90 (ko; u=1,--+,m), 
and, since the determinant in (4) does not vanish, 
x, +0 (kR— o;p=1,---,m). 


We see that the modified procedure is indeed convergent. 

10. In comparing the rate of convergence of the original and the modified 
procedure it is better to change our notations in the following way. If we start 
with a vector £) = Xo and apply the complete -cycle of single steps correspond- 
ing to N;, = 1, ---, , the obtained vector will be denoted by X;, and the vectors 
obtained in repeating each time the complete n-cycle will be denoted by X2, 
Xe -++, Xn = by °° 

In the same way, in the modified cyclic procedure we obtain, starting from a 
vector {9 = Zo and applying each time the whole (n + 1)-cycle corresponding to 


Ni; = 0,1, ---, , the sequence of vectors Z;, Z2, ---, Zi = fnk, ***- 
11. We decompose A in the following way: 
(15) A=L+D41"*, 
where D is the diagonal matrix with the elements a1, ---,@nn, while in L all 


elements to the right and on the main diagonal and in L* all elements to the left 
and on the main diagonal vanish. Then the rate of convergence of the usual cyclic 
single step iteration depends on the maximum modulus Ay of the roots of the 
equation 


(16) IA(L + D) + L*| =0, 
and we have, if Ay > 0: 


(17) Xi = O(ntk™) (k- &), 


while the starting vector X» can be chosen so that X;Aw~* does not tend to 0 as 
k — o. The proof of this is quite similar to the following proof of the correspond- 
ing results for the modified single step iteration; see (25). If Aw = 0, then X; 
vanishes identically, and the solution is obtained at the most in m steps. 

12. We will now characterize in a similar way the rate of convergence of the 
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modified procedure. We have for the matrix A of the system (10) the decomposi- 
tion corresponding to (15): 


(18) 4-(% ww) aL+D+L 


and obtain between Zp» and Z,, as in the theory of the usual cyclic one step 
iteration, the relations 


(19) (L + D)Z, + 1*Z, = 0, 
Z,= — (L+ D)2t*Z,. 

13. As has been mentioned above the result of this operation is not changed if 
Zo is replaced by the corresponding reduced vector 25 = (0, x1, ---, x,). 
Before we go on from Z;, we replace therefore Z; again by the corresponding 
reduced vector 2; = (0, x:™, ---,x,™). For this purpose we apply the trans- 
formation x, = Z, — Zo (u = 1, ---, ) which is equivalent to left multiplication 
by the (w + 1) X (m + 1)-matrix 


(20) m=(% 2): 


where the first row consists of 0’s and the first column, with the exception of the 
0 on the top, of —1’s. 
We have then finally, putting 


(21) Q=- NAL + D)=2*, 
(22) 2. = OckZo (k = 1,2, ---). 
14. If Qo vanishes identically, we have at once Z; = --- = Z, = --- = 0,and 


the solution of (1) is attained at the first step. We can and will therefore assume 
that Qo + 0. We use then the following result due to Werner Gautschi [5]. 
If for any matrix C = (c,,) we define as its “norm” 


N(C) = vy | Cus |?, 


then if C is a square matrix of order m for which the greatest modulus of a funda- 
mental root is A, we have 


(23) N(C*) = O(A*ke) (kR— &), 


where p is the greatest multiplicity of a fundamental root of C with the modulus A. 
As a matter of fact (23) is not the best possible result, cf. [10], p. 5, Satz V. But 
Werner Gautschi’s result is completely sufficient for our purpose. 

15. If we apply this to the singular matrix (21) and denote the maximal 
modulus of a fundamental root of Qo by Ag, then p does not exceed n — 1 if 
Xe > 0, as will follow later from (31). We have therefore 


(24) N(Qo*) = O(\e*k™*)  (k- @). 
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Further it follows from (22), in applying the Cauchy-Schwarz inequality, 


|2e| < N(Qo®)|Zol, 
and we obtain therefore 


(25) 2. = O(rctk") (k- @). 


16. On the other hand it is easy to show that for a conveniently chosen starting 
vector Zy the expression Z;Ag~* does not tend to zero. Indeed, if n is an eigenvector 
of Qo corresponding to a fundamental root \ with |A| = Ae, we have 


An’ = Qon’ 
and iterating 
A‘n’ = Qokn’. 


But, since the first row in Qo consists of zeros, the vector 7 is a reduced one and 
can be taken as 29. Then we have 


: r \F 
2.= 20, 2Zre* = (>) 20, 
Xe 
and this does not tend to zero as k > ~. If Xg = 0, then Z, vanishes identically. 


17. We shall now transform the fundamental equation of Qo, and we introduce 
for this purpose the matrix 


- m= (4. 2), 


where the first row and the first column consist respectively of 0’s and —1’s with 
the exception of the first element e. NV, corresponds to the transformation 


(27) yo = Lo, y = Z,— Zo (vy = 1, +++, m) 


and for ¢ — 0 goes into No. Since the inverse of (27) is for « ¥ 0 


1 1 
Zo = = Io Z, = Y += yo (v = 1,---,m), 
we have 


(28) No = a Me 


where the first column consists of «~', while the first row with the exception of 
the first element contains only 0’s. 
The fundamental equation of Qo can be written in the form 


(29) lim |AE + N.(L + D)=7L*| = 0. 


18. On the other hand we have identically, since | V.| = e, 


(30) (J2L+D))(IAE+ NAL + D)7L*|) = ela(l+ Dn + LI, 














of 
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and we obtain the fundamental equation of Q» in taking the limit for « ~ 0 on 
the right in (30). 
Now we have 


we € doo 0 
(+ D)N& = (Se L+ a 


where Gno + Gai + -** + Gann = 0 by (9); and we have therefore identically 


doo aor 


4, A(L+D)+L* 


Ag is therefore the maximum modulus of the roots of the equation of degree n 


e[A(L + BD) + 2*|=r 








-| oo do» 
Ir A(L + D) + L* 


One root of this equation is zero. 
19. In specializing for » = 2 we obtain in particular, if we put a,,; = a, 
22 = G2, G12 = G2; = o and assume o + 0: 


(31) G(A) = 0. 














@,+a2+20 —-adi-—0 —G@2-<¢6 
G2(A) = —o a, a 
0 Xo Ade 
= (a; + a2 + 20)aea;A? — od(a1 + oc) (a2 + 0), 
(32) oe |o(a; + o)(a2 + @)| 


@102(a; + a2 + 20) ’ 


while the equation (16) for Aw reduces to \*a,a2 — Ao? = 0, and gives 





(33) Ay = 
From (32) and (33) we have 


Ae |aite|laz+o| 


(34) Aw [o| (a, + a2 + 2c) 





20. If we square this, subtract 1 and multiply by the square of the denomi- 
nator we obtain 


(a: + o) (a2 + ¢) — o(a; + a2) — 207 [ (a1 + 2) (G2 + &) + o (a; + G2) + 207] 
= (@:@2 — 0”) (aya2 + 2(a; + a2)o + 3o%). 


Since the first factor is positive, we see that Ag = Aw according as 
(35) o(c) = 30? + 2(a, + as)o + aya, S 0. 


Here ¢ ¥ 0 is subject only to the condition |o| < Vaya2. We have obviously 





o(—Vayar) = 4Va,a, (aw, ~% + *), 
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4 % Be : a\a 
Since this is <0, we see that ¢(—c) has two positive roots oi, o2 with o,02 = = 4 


and 
0 < a; < Vaj2 K oz. 


It follows that for = 2 we have Xg < Aw, if o1 < —o < Va,a2, and Ag > Ay, if 
~ Vaya2 < —o@ <;. The modified procedure is in the first case faster and in the 
second slower than the original one. 

21. To prove the corresponding result for m > 2 consider the matrix A of the 
quadratic form 


(36) K(é) = ayxy? + 2oxyx2 + aoxe? + D x,?. 


u=3 
In the corresponding determinant (31) for G(A) the elements in the first column are 
Gyo + Oui + °°? + Oy = — (Gy, u41 + -+* +b Ops) 


and vanish therefore for u > 2. The same is true for the elements \a,, to the left 
of the main diagonal with »> 2 and v < yu, while the elements on the diagonal, 
Ady, (u > 2), become respectively \. We obtain therefore 


G(A) = A**G2(A) 


so that Ag is given in this case by (32). 

22. In the same way it follows from (16) that Ay is again given by (33). We 
can have in this case, according to the chosen values of o, either X\g > Aw or 
Age < An. - 

It may finally be remarked that the value of \g is not changed, if the (m + 1)-st 
equation in (10) and the corresponding new variable 2» are not put at the beginning 
but are interpolated between two indices u, w + 1 or even put at the end. Indeed 
this amounts to the old process applied to a transformation of & by a finite 
sequence of single step iterations, but then £ is carried over into the general 
reduced vector, and the invariance of \¢ follows then from the characterization 
of \g contained in the developments of numbers 15 and 16. 

23. We consider now Seidel’s relaxation procedure (8). Then in the 
modified procedure we obtain the speeding up for an index k for which we have 


(37) (ro™)?/Ao > (r4)?/A, (u #9), 
that is to say, by (12): 


(38) [E16 |/Aat> [74 /A,M(ux60). 


y= 


We will now show that this inequality is impossible, if we have 


(39) Aé > DA, — min A,!. 
v= l<gucn 
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Indeed, if we put 
|r, | = p,A,}, = 1, ---,n, 


we have from (38) and (39), since one of the 7, and therefore one of the , vanish: 


[Xo 75™|/Act < X p,(A,/Ao)t < Act(— min A, + > A,}) max p, 
y=1 v=1 l<egn v=1 
€ max p, = max (|7,|A,-). 
l<eqn 


Therefore, in order that (37) be possible at all we must have 


(40) V4.< 5 V4.- min VA, 


» 
u=l l<gucgn 


24. In order to discuss the situation under the condition (40) we consider the 
r,™ as stochastic variables and discuss the probability for (37) under suitable 
assumptions on the distribution of the r,™. We put 


(41) 6,= VA,>0 (vy =0,1,---,m) 


and denote the sequence of the vectors p(k = 0,1, ---) by S. In the following 
formulae the index of p, that is the upper index of r,™, will be dropped whenever 
it is possible without danger of misunderstanding. 

Then our problem can be described as the problem of computing the probability 


(42) P| 


defined in the usual way as the limit of the relative density, and ascertaining 
whether (42) is positive. As the relations in the brackets of (42) are homogeneous 
in the 7,, each of the vectors p e S can and is from now on assumed to be normed 
in such a way that we have 

Ir, 





', 
> max | l. pes| 


l<qrqgn sd 


¥ 


v=] 








(43) max —— = 1. 

l<grgn B, 

25. Denote the N-th section (¢;, ---, pw) of S by S™ and put 
(44) #0) = P[ En <o, ve sm]. 
v=] 
This is the “finite” probability in the classical sense, and we put then 
(45) @(c) = lim 6 (¢), 
Noo 


if this limit exists. 
For any k, k = 1, ---,, we denote by S, the sequence of the vectors p,, 
\ < N, for which k is leading index corresponding to p,_; and therefore 7. = 0. 
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Further we write S, for the complete partial sequence of S corresponding to 
these vectors. 

We denote further by Z,(NV) the number of vectors p, in S™ whose leading 
index is » and put 


(46) }.%)(¢) = P | Sr <o; pe sn]. 


v=] 


Denote now by S;,, the partial sequence of S containing the vectors p, such 
that yu is leading index for p, and & is leading index for p,_;, and therefore by (8) 
and (43) r? = 0, |r® | = B,, and by SY? the sequence of the vectors p, from 
Sk,» with A < N. Put then 


(47) 22(0) = P[ En < ones], 


and denote by Hi.,,(V) the probability that the leading index u in S™ follows 
the leading index k. 
26. We define further 


1 

(48) HM, = lim yy Le), 

(49) Hi, = lim Hi, (N), 
No 

(50) &,(c) = lim 69 (ce), 
No 

(51) @,,,(¢) = lim $f? (c) 
Noo 


if the limits (48)—(51) exist. We have then obviously by elementary probabilities 


— + o™,.) OY 
(52) (0) = LEP) Ta 
(53) oP (o) = L Of) (¢)M..(N). 


We see that, if the limits (48), (49) and (51) exist, &,(¢) and 6(c) also exist and 
we have 


(54) 5(c) = ¥ 4(o) Hs, 
(55) &(c) = > , (0) Ak, p. 


27. We define now m vector-fields F, consisting respectively of all vectors p 
normed according to (43), for which the component r; vanishes. Each of the fields 


~ 


WAT... eee 
x ae 
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F, can again be decomposed into m — 1 vector-fields F;,,(% ~ «) consisting re- 
spectively of all vectors p normed according to (43) with r, = 0, |r,|= B,. A 
vector p can of course belong to several of these fields. 

In the fields F;, F;.,, we will now define in a suitable way the following geometric 
probabilities: 


(56) nite) = P| Dn < eich], 
y=] 

(57) wr» = Pl |r| =B,; pe Fel, 

(58) $*,.(¢) = P [= r Qojpe Fas]. 


Of course, to define these geometric probabilities we must use a convenient model, 
which presents itself here in a very natural way. If the components r,™ of the 


residual vector p, are computed up to a certain number of decimal places, they 
a) 

are approximated by certain fractions a , where m is a power of 10. It is then 

m 


natural to assume that, if the subscripts k and u are fixed, the integers g,™ vary 
uniformly within their limits. We consider therefore m + 1 integers 


m, mM, (vy = 1, «++, m) 


tending to infinity in such a way that we have 


pa a. Soe Y 
m 
We define then 
(59) &*(o) = lim P| Ea so ilies Rear a 0], 
mo v=1 l<v<n My, 
saith sein tei NA ce ha ht 
(60) wr,» = lim P | qu| =m,; Max =i,qg=0], 
mo l<qgrgn m, 
(61) #4 = lim P| Ea < om; las] = mu ae = 0, max /@l - 4]. 
m->0 v=1 l<grcgn m, 


28. It is easy to compute the w,,,. Assume k = 1. Then the probability o,, 
in (60) is obtained in considering the right parallelepiped |x| < m, (A = 2, ---, #), 
as the quotient of the area of the two faces |x,| = m, by the complete area of all 





faces. In putting temporarily m2 --- m, = M we obtain 
_ 2M /m, _ m/m, 
22 M/m YL m/m 
A=2 A=2 
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and this tends to 
1/8, 


> 1/8, a 1/8; 


v=l 


In replacing here 6, by 6; we obtain finally 








(62) Ok» = = /B» (k ¥ u). 
i 1/8, a 1/6; 


v=l1 


We write this differently in introducing the expressions 








1 
(63) oo. = a. 
X 1/6, 
v=1 
Then we have obviously 
(64) > a = ‘, 
A=1 
and (62) becomes 
Cn _ oy 
n ‘ats 1 a 3 
Lo — ” 
v=] 
If we now define w, as 0, we obtain 
on 
k , 
(65) a {i = = ( x ») 
0 (k = p). 


The values ©*;,,(¢) and &,*(c) have already been computed in another com- 
munication [8], and we will give *;,,(c) later on. In particular we have 
(66) ,* (c) = Wk, p®* x, u(o). 
nak 


29. We make now two fundamental assumptions about the sequence S under 
consideration : 


(67) Hypothesis A. Gup(0) = O*np(0) (kX ui kyw =1,---)m) 
(68) Hypothesis B. Hi,.(N) > we,n (N+ ©,k ¥u;k,u =1,---,n). 
It follows then from (55) that we have 


(69) 4, (co) = &*(c). 


In order to obtain the value of &(c) we will now prove the existence of the 
limit (48). 
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30. We have obviously 


(70) > Li(N) = N +1, 
k=l 
L,(N) _ a Ai, (N)Li(N _ 1). 
ime 
If we replace here each L;(N — 1) by L;(N), the modulus of the error does not 
exceed >> Hi, = 1, and we obtain 
ied 

DY Aiw(N)Li(N) = L,(N) + 4(N), |e(N)| < 1. 

ime 
Divide this by N and put 
L,(N) 


(71 
(71) N 


= h,(N). 
Then we obtain 
ms 1 
(72) Y Hew(N)b(N) — (N) = & (NY), «(N) = o(+) (u = 1,---,m). 
in 
In the same way we obtain from (70) 
“ 1 
(73) DL m&(N) =1+—- 
k=l N 


If we replace in the equations (72) and (73) the h,(N) by x,, the equations become 
for N — « by (68) and (65) 





= x, 
(74) + + = % (u = 1,---,m), 
kat l — o% 
ku 
(75) > x, = 1. 
k=1 
We consider first the system (74). If we divide by o, and add = on both 
— 


sides, we obtain 


Xk Xu 





= Xn + Xn a 
pai l — o% on 1—o, o,(1 — o,) 





Denoting the expression on the left by y we get 
(76) X, = you(1 — oy). 


31. On the other hand, if we put the value (76) for y = 1 into the system (74), 
the system is satisfied; we see that the rank of (74) is exactly m — 1. We can 
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therefore choose n — 1 of the equations (74) in such a way that they form, taken 
together with (75), a system of the rank m. Therefore the corresponding n — 1 
equations (72) form, taken together with the equation (73), a linear system which 
remains regular in the limit N — . Therefore its solutions 4,(N) tend to the 
solution of the system (74), (75); and this solution is obtained from (76), if y 
is chosen so that (75) is fulfilled. But then we have 


I 


<%= 1 =r(1- Eo), 
p=l1 w= 


(77) Li.(N) _ on(1 — ox) 


1 — >a,’ 


y= 


If we introduce here the values (63) of the o,, we obtain 
Be By — By) 
208% 


vr 





(78) i, = 


It follows now from (54), (55), (68), (65): 


(79) @{e) = x A = ton, ,B*s,»(0). 
nat 


32. An explicit formula for &*;,,(¢) can be written simply by using the expres- 
sion for a function F(¢) which was defined and computed in a previous paper [8 ]. 
In replacing the integer ” used in that paper by m — 2 (for m > 3) we define 
F(c) as the probability 


(80) Fo) = Plé&i +f +--+ tone Ko; |&|Ka (vw =1,---,m — 2)]. 


The a, are » — 2 positive numbers. Then we have for F(c) (see formula (5) 
in [8 ]): 
1 1 —s 


(81) F(c) = : II (1 — S)(a@ + o)5’, 


2"? (m — 2)! ar +++ Ga-2 et 
@=ai+ +++ + an2. 





Here S is the operator defined by 
S'f(a) = f(a — n). 
The computed expression consists of 2"-? terms of the form 


+ (a — 2m, — 2m, — --- +0)”. 


The subscript + signifies that, if the expression within the parentheses is negative, 
the whole expression has to be replaced by 0, while otherwise the subscript can 
be dropped. 

In particular we have pointed out in [8], p. 6, that F(c) is strictly monotonically 
increasing with o for — a < o < a and grows from F(—a) = 0 to F(a) = 1. 

















ON GAUSS’ SPEEDING UP DEVICE 131 


We see that we have 
(82) 1>F(e)>90 (je|<a). 


F(c) is further a continuous function of ¢. We obtain therefore the same value 
of the probability, if we replace in (80) the condition <o¢ by the condition <ce. 
33. &*, .(c) is by definition (58) the arithmetical mean of the two probabilities 


P| os éo~hsleie e,|.P| PET AT |. 
a, ot, 


We have therefore 
(83) $;, .(c) _ $L Fin (o — By) + Filo + B,) J. 


The function F;,, is here obtained from (81), if the a;, ---,a,—2 are identified 
with the » — 2 of the numbers 8, which remain after deleting 8, and 6. Intro- 
ducing (83) into (79) we obtain finally the complete expression of @(c). 

If we return now to the problem of computing the probability (42) and norm 
here the 7, according to (43), we see that we have to compute 


P| > Bo, max Ir. = 1, min |r,|= 0}. 


l<grvqgn ’ lgrqn 
This is the sum of the two probabilities 


n 
2X 


v=1 

















P| Er. < — fo; max [7 = 1, min ini= ol, 
yan] l<grgn B, l<grgn 
or ans eae SS 
P| > r, > Bo; max = 1, min |r,|= 0]. 
v= lgvrgn PP l<grgn 


But these two probabilities have the same value, and the first of them is obviously 
@(—£o). We obtain therefore for the probability (42) the expression 


(84) 26(—£o) = 2 Ay u wr, wl Fr,u(— Bo — By) + Fi,n(— Bo + By) J. 
unk 


34. We will now prove that (84) is positive under the condition (40), that is 


(85) Bo < £8, — min ,. 
v=] l<grgn 


Observe that all terms in (84) are non-negative. It is therefore sufficient to prove 

that under the condition (85) at least one of the expressions F;,,(— 80 + 8,) is 

positive, that is to say, by (82), that for convenient choice of k and yu either 
n 

Bo + 8, or |Bo — B,| is less than >> 8, — 8, — Bx. Assume now that the 8, are 


vy=1 


ordered increasingly 
Bi gf B2Q +++ LBs 
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and put 
B= 686i + --- + Bs. 


Then I say that for k = 1, 4» = 2 we always have 


(86) |Bo — B2| < B — Bi — Ba, 
that is to say that 
F1,2(— Bo + B2) > 0. 


To prove (86) consider the two cases By > Bz and Bo < fe. In the first case (86) 
reduces to 


Bo — B2 <8 — Bi — Bx, Bo < 8B — Bi, 
and this is exactly (85). In the second case (86) reduces to 


B2 — Bo <8 — Bi — B2, Bo > — (8B — Bi — 282); 
but here the expression on the right is not positive, as 


(Bs — 82) +Bt+- + + Bn 2 0. 


We have now proved that under the condition (40) and the assumptions of 
Nr. 29 the probability for the speeding up of Seidel’s relaxation procedure in using 
Gauss’ device is positive. 

The author gratefully acknowledges discussions with T. S. Motzkin and 
O. Taussky-Todd. 
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Extension of Lindow’s Tables for Numerical 
Differentiation Using Newton-Stirling and 
Newton-Bessel Differences 


By Herbert E. Salzer and Genevieve M. Kimbro 


Lindow [1] gives several short tables of coefficients for obtaining the first 
and second derivative at intermediate points by differentiation of either the 
Newton-Stirling or Newton-Bessel interpolation formula [2]. The former involves 
central differences uio?"*! and 69" which are on the line with f(x»), and the latter 


involves central differences 6,;2"** and yé,;" which are on the line with f (xe + =) - 


h being the tabular interval. Lindow’s tables of Newton-Stirling coefficients are 
intended to interpolate for f’ (xo + ph) and f’’ (xo + ph) for either 0 < p < .25 or 
for .75 < p < 1.00, in the latter case by choosing a new x» = original x» +h 
(and corresponding central differences opposite the next entry) and a new p 
where now — .25 < new  < 0. His tables of Newton-Bessel coefficients are in- 
tended to interpolate for f’ (xo + ph) and f” (xo + ph) for .25 < p < .75, or in the 
new variable p; given by p = 4 + pi, for — .25 < 1 < .25. The reason for 
employing both Newton-Stirling and Newton-Bessel formulas (i.e. both tabular 
and mid-point central differences) is that by proper choice of p or p; the argument 
is never more than one-fourth of a tabular interval away from the central differ- 
ences employed, which should yield high accuracy. 
The formulas for f’(x) and f” (x) are of the form: 


(1) f’(xo+ ph) = ; {wb Pb + YE (Asmar (P) ude + Aomso(p)be"*) | + Ri, n 
m=1 


obtained by single differentiation of the Newton-Stirling interpolation formula, 
to be used for + ranging from —.25 to .25. 


(2) f'(xo + ph) =f (x. + ; a pik) 


1 nn 
=3% {a + pydy + DO (Bomgi(h1)dpe" + Bomyo(P1) mdz") | + Rens 


m=1 





obtained by single differentiation of the Newton-Bessel interpolation formula, to 
be used for p = 4 + p1 ranging from .25 to .75, or +p; ranging from —.25 to .25. 


(3) f" (xo + ph) 


. a 4 pabé +E (Comsa(p)i Comsa(P)ade*)| + Raw 


m=1 


obtained by double differentiation of the Newton-Stirling interpolation formula, 
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to be used for +p ranging from —.25 to .25. 
h 
4) f" G0 + ph) =f” (xo +5 + ph) 


1 n 
= ie {ay + pid + YD (Domy2(1) mix"? + Domss(p1)5;"**) | + Ra ay 


m=1 


obtained by double differentiation of the Newton-Bessel interpolation formula, 
to be used for p = § + pi ranging from .25 to .75, or +; ranging from —.25 to .25. 

In the mth pair of terms only, the right hand term within the parentheses does 
not appear if the interpolation series (1), (2), (3), or (4) is not taken beyond 
pier"*, 62"41, 59?**? or wd,?"** respectively. The formulas for the remainder terms 
Ria, t = 1, 2, 3, or 4 depend also upon the parity of the order of the last difference 
retained in (1), (2), (3), or (4). To find R;, explicitly one may differentiate the 
formulas for the remainder terms in the Newton-Stirling or Newton-Bessel 
formula ([3], p. 100, 102), making use of the properties of divided differences 
under repeated differentiation ([3], p. 66-67). For (1)—(4) there will be altogether 
eight different formulas for R;,,, having either two, three, four or six terms, the 
use of which is quite troublesome and time-consuming, since each term involves 
the evaluation of a polynomial of higher degree and the estimation of a higher 
order derivative. For most practical problems one can almost always bypass the 
work in using the theoretically exact expression for R;,,, by simply observing 
the last term retained in any of (1)—(4) in conjunction with the rate at which the 
terms are falling off in magnitude. 

Following are the explicit expressions for the coefficients A,(p), B,(p1), C,(p) 
and D,(p), as (even or odd) polynomials in p or ;, for r = 3(1)10 and s = 4(1)10: 





























A;(p) = ee L. A,(p) = eae As(p) = 5p! ot + 4 
esate te Mere ue +” , gc apt a Atte = 36 
Lins Eo Se 16 — 2460p" + S76 

Aw(p) = 5p® — 120p’ + Sie ene + 576p 

Bi(p,) = PE, Bu(p,) = “PSP 

By(p,) = 20P# ~ ee $9 5 See ss + 2591 
B,(p,) = Atbes* = 2800pi* + 3108p," — 225 





3 22560 a 
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Ba(p,) = 64p17 — 1008p,° + 3948p,* — 3229p, 
oP) = 3 22560 : 
ited 2304p,* — 37632p,° + 1 57920p,4 — 1 54992p,? + 11025 
(Pi) = 928 97280 ; 
ise 1280p,° — 42240p,7 + 4 21344p,5 — 13 82480p,* + 10 57221p, 
0(Pi) = 4644 86400 
_ 6p -1 _ 2p* — 3p 
Cs(p) a 12 ’ Cs(p) — 12 ’ 
_ 15p* — 30p? + 4 _ 3p5 — 20p* + 21p 
_ 14p* — 105p* + 147p* — 18 
_ 6p? — 105p* + 455p* — 410p 
Gis 15p* — 280p* + 1365p — 1640p? + 192 
o(P) = 6 04800 
126° — 5 4p? — 3 
D4(p1) -_-—., Ds(1) ced 
_ 240i — 840p;* + 259 _ 336p:5 — 1400p + 777p, 
Ds(p1) = 5760 » Dilb:) = 40320 ; 
448p,° — 5040p; + 11844p,2 — 3229 
Ds(b1) = 3 22560 
A _ 5176p: — 7056p:5 + 19740p;* — 9687p, 
o(P1) = 29 03040 ¥ 
- _ 11520p:* — 2 95680p:° + 21 06720p;* — 41 47440p,* + 10 57221 
10(P1) = 4644 86400 ; 


The purpose of these present tables is to extend Lindow’s tables which go 
only as far as the coefficients of the 6th differences in formulas (1)—(4) above, 
and give only 5D up to the 4th difference coefficients and 4D for the 5th and 6th 
difference coefficients. These present tables give the coefficients A,(p), B,(p;), 
C.(p) and D,(p1) at the same interval of .01 and the same range of p or p; from 0 
to .25, as occur in Lindow, but for every difference as far as the 10th difference 
inclusive, and to 10 significant figures. This represents a considerable extension 
of Lindow’s original tables and should be useful in many calculations of the first 
and second derivatives at intermediate points, which arise in numerical differenti- 
ation work and in the solution of first or second order differential equations, where 
Lindow’s tables are entirely inadequate. Anyone who has performed numerical 
differentiation, especially for the second derivative, is aware of the great loss in 
significant figures due to the power of h in the denominator, as well as the sub- 
traction of nearly equal terms in the numerator. For functions that are not deter- 








136 EXTENSION OF LINDOW’S TABLES FOR NUMERICAL DIFFERENTIATION 


mined by either measurement, observation, experiment or approximation, but 
which are mathematically defined so as to be computable to any degree of pre- 
cision, the only limitation to the accuracy in numerical differentiation is due to 
the truncating error in using a finite number of terms of the formulas (1)—(4) 
and the computing error due to the carrying of a fixed number of decimal places 
or significant figures in the computation. Thus Lindow’s original tables severely 
limit the accuracy attainable for mathematically defined functions because of 
those two mentioned reasons. These present 10S tables, as far as the 10th differ- 
ence, are intended primarily to reduce considerably both truncation and com- 
putational errors. 

The only other conveniently available tables for numerical differentiation at 
intermediate points, employing central-type differences, appear to be those of 
Davis [4] who gives, at intervals of .01, the Newton-Stirling and Newton-Bessel 
coefficients, but only as far as the 5th difference and only to 5D, and the Everett 
[2] coefficients as far as the 6th central difference and to 10D. However, all of 
Davis’ coefficients are for the first derivative only. Davis’ 10D Everett coefficients, 
giving 8-point accuracy, might suffice for most problems requiring just the first 
derivative. But Lindow must be extended anyhow to take care of the equally 
important second derivative, and thus we might as well be uniform in procedure 
and accuracy and use Lindow’s arrangement for finding also the first derivative. 

The present calculation was done originally using only a desk calculator, by 
exact computation of the numerators in the coefficients A,(p), B,(p1), C,(p) and 
D,(p1) in (1)-(4), and rounding only after division by the denominators. Thus 
all tabular entries should be correct to within a half-unit in the last (tenth) 
significant figure. All entries on the preliminary manuscript were rechecked by 
Norman Levine, employing the IBM 704, using double-precision (floating point) 
arithmetic. An additional differencing check was perforrhed by hand upon the 
entries on the preliminary manuscript. Also a functional check was performed 
upon every entry in the typewritten final manuscript by computing the first and 
second derivatives of (1 + x) for x = 0(.01).25 and x = .50(.01).75. 


Convair-Astronautics 
San Diego 12, California 
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2. E. T. WHITTAKER & G. Rosinson, Calculus of Observations, 4th ed., Blackie & Son, London, 
mee p. 38-41. 
F. B. HILDEBRAND, Introduction to Numerical Analysis, McGraw-Hill, New York, 1956, 
p. 66-63, 99-103. 
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TECHNICAL NOTES AND SHORT PAPERS 


Table of Integers Not Exceeding 10 00000 That Are Not Expressible 
as the Sum of Four Tetrahedral Numbers 


By Herbert E. Salzer and Norman Levine 


For the past 21 years one of the authors has been concerned with empirical 
theorems expressing certain classes of positive integers as the sum of four tetra- 
hedral numbers, or “‘tetrahedrals’’ T, = n(n + 1)(m + 2)/6, n > 0, the results 
being contained in short notes and abstracts, as well as unpublished tables 
[1}+[9]. We use 2, to denote a sum of p T,’s. Among the move interesting past 
findings were: every square < 10 00000 a 2, [7], [8], every T, < 10 00000 a =, 
other than the trivial decomposition 7, = T, +0+0+0 [3], [4], every 
multiple of 5 < 1 00000 a 2, [9] (see also [11]), and verification of Pollock’s 
conjecture that every integer is a 2; for the first 20000 integers [5], [10]. Inci- 
dentally all investigations in [1 ]}-[9] were done by hand, employing at the most 
a desk calculator. 

The exceptional numbers, which by definition are those not expressible as 2,, 
were tabulated previously only up to 2000. But even that far interesting features 
turned up, such as 1314 being the only exceptional number ending in 4, and the 
very few ending in 1 or 9 [4]. Then for numbers ending in 6 and < 20000, 6186 
turned out to be the only exceptional one [6 ]. Also there was no striking difference 
between the density of exceptional numbers in the first thousand and in the 
second thousand brackets, decreasing from around 44% to 3%, so that it was 
interesting to speculate upon the approximate density in the neighborhood of say 
10 00000. Then finally, the verification of the conjectures that every integer is 
a =; and that every integer m = 5r is a X, for the first 20000 cases, while every 
10r + 6,7 = 0,1, 2, ---, 617, is a 24 until 6186, posed the question as to whether 
an exception might occur to either of the former two empirical theorems even 
after verification in those first 20000 cases. Thus it was felt that tabulation of the 
exceptional numbers < 10 00000 would afford a much clearer picture as to their 
distribution and density, as well as stronger evidence for the truth of Pollock’s 
conjecture and the author’s (and Richmond’s [11]) conjecture that every 
m = Srisa X4. 

This table of exceptional numbers < 10 43999 presents a great surprise in its 
picture of their distribution which is entirely different from that envisioned from 
those < 2000 [4]. Most strikingly unexpected is the decrease in the density of 
exceptional numbers from several percent in the neighborhood of 2000 to what 
appears to be practically zero near 10 00000. The scarcity of exceptional numbers 
in the higher ranges is in accordance with Hua’s result that ‘almost all’’ positive 
integers are expressible as a 2, [15], [16]. But even more, this table shows that 
the likelihood of a given number m being exceptional falls off so rapidly with 
increasing m that it appears to be a plausible conjecture that there might be some 
my sufficiently large such that every m > mp is a 24 (a conjecture unlikely to sug- 
-gest itself from the exceptional numbers among only the first few thousand). 


Received 14 January 1957. 
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TABLE OF EXCEPTIONAL NUMBERS < 10 00000 
17 1227 3183 9772 29157 € 
27 1233 3218 9973 29487 1 
33 1243 3263 10397 29938 
52 1314 3463 10467 30298 
73 1382 3512 10532 31973 
82 1402 3887 10633 33183 
83 1468 4003 10852 36262 
103 1478 4307 11237 36913 
107 1513 4317 11302 37798 
137 1523 4563 11737 38453 
153 1578 4832 11962 38707 
162 1612 4923 12247 38807 
217 1622 5013 12547 39693 
219 1658 5142 12722 39913 
227 1678 5238 12777 41278 
237 1693 5283 12843 41322 
247 1731 5483 12858 41433 
258 1738 5508 13127 44833 
268 1742 5538 13393 47627 
271 1758 5563 13822 48043 
282 1767 5618 14492 56467 
283 1803 5647 15122 56842 
302 1858 5707 15483 58613 
303 1907 6022 15867 59077 
313 1923 6057 16097 62158 
358 1933 6067 16538 64752 
383 2037 6186 16637 65253 
432 2053 6213 16742 - 65567 
437 2172 6263 17253 71157 
443 2198 6343 17683 74687 
447 2217 6462 17813 78003 
502 2218 6863 17893 78787 
548 2251 7067 18573 83603 
557 2253 7278 18782 84023 
558 2327 7377 19168 85993 
647 2372 7387 19277 91128 
662 2382 7423 20918 1 06277 
667 2417 7497 21523 1 13062 
709 2437 7542 22618 1 34038 
713 2457 7662 22657 1 48437 
718 2537 7793 23677 3 43867 
722 2538 7873 24237 
842 2578 8223 24317 
863 2687 8307 24338 
898 2818 8322 25447 
953 2858 8973 25723 
1007 2898 9063 26007 
1117 2973 9488 27858 
1118 3138 9687 28617 
1153 3142 9753 28847 
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A quick glance over this table is sufficient to verify Pollock’s conjecture for 
each of the first million integers. Furthermore, from the absence of any excep- 
tional numbers in the range 3 43867 to 10 00000 and from an upper bound to 
the magnitude of the first differences of tetrahedral numbers < 2500 00000, it is 
easily shown that for any number m between 10 00000 and 2500 00000, it is al- 
ways possible to find a T, such that 3 43867 < m — T, < 10 00000. Thus every 
m < 2500 00000 is a 2s. (This fact and the extreme rarity of large exceptional 
numbers makes the empirical theorem of Pollock that every integer is an 2; as 
overwhelmingly certain as one can be short of actual proof.) This table verifies 
that every m = 5r is a 2, for the first 2 00000 values of r, making it also a very 
plausible empirical theorem. In the class of the first 1 00000 values of both 
m = 10r + 4 and m = 10r + 6, 1314 and 6186 respectively are the sole excep- 
tional numbers. Among the first 1 00000 values of m = 10r + 1, only 271,1731 
and 2251 are exceptional. Among the first 1 00000 values of m = 10r + 9, only 219 
and 709 are exceptional. All exceptional numbers m, where 6186 < m < 10 00000, 
are here seen to be of the form 10r + 2, 107 + 3, 107 + 7 or 107 + 8, and that 
appears to be another plausible conjecture for every exceptional m > 6186. 

The method of computation was to obtain (Pe) +{**),('47) 
for every 1, j, k, | > 0 until one found for every m < 10 43999 either a represen- 
tation as a 2, or that that m was exceptional. The calculation was begun upon the 
Univac Scientific Computer (ERA 1103) at the Convair Digital Computing 
Laboratory, the initial part re-run, and then those results were checked and 
continued upon the IBM 704 Digital Computer. At the start upon the ERA 1103, 
500 words of high speed storage with 36 binary bits in each word permitted the 
investigation of 18000 numbers at a time. Then the IBM 704 had at its disposal 
3000 words of high speed memory, each of the 36 binary bits in a word represent- 
ing a number, so that 1 08000 numbers could be investigated at one time. The 
1 08000 binary bits were filled with 1’s at the start, and a 0 introduced into the 
binary position of the word which represented a non-exceptional m. After all 
possible cornbinations of i, 7, k, } had been exhausted, each of the 3000 words was 
searched for binary bits that remained 1 and the exceptional numbers m corre- 
sponding to those bits were printed out. In choosing combinations of i, j, k, I, 
repetitions due to symmetry were avoided, as well as combinations yielding an 
m that was either too large or too small for the group of 1 08000 numbers under 
consideration. 

Those interested in actual mathematical proofs (which appear to be rather 
involved) may consult Dickson [12], [13] for earlier work, and Watson [14] 
for the sharpest results to date. It is rather amazing that the proved p in m = 2, 
for every m is p = 8, and no better than = 8 for arbitrarily large m, while the 
actual p (according to the evidence in this table) may be only 5 for every m, 
and 4 for sufficiently large m, less by 3 and 4 respectively. Considering the diffi- 
culty of the existing proof for p = 8 [14], one may well wonder, should p = 5 or 
pb = 4 be the truly minimum values, for every m, and m sufficiently large, respec- 
tively, how long the world must wait and how difficult and sharp the mathematical 








144 GROUPS OF PRIMES HAVING MAXIMUM DENSITY 


tools must be, until the desired proofs would be found. 
References [15 ] and [16 ] were called to the author’s attention by K. A. Hirsch. 
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GROUPS OF PRIMES HAVING MAXIMUM DENSITY 
By John Leech 


The following lists give groups of six or more primes which minimize the 
difference between first and last, the lists being complete for the range 50 to 
100 00000. Four numbers out of nine can be prime, such as 191, 193, 197, 199. 
There are 897 such groups of four in the range. Five numbers out of thirteen can 
be prime; there are 318 such groups in the range. Six numbers out of seventeen 
can be prime, such as 97, 101, 103, 107, 109, 113; there are seventeen such groups 
in the range, centered on: 


Received February 7, a Due to misfiling in the MTAC office, this paper is appearing later 
than was scheduled ; see MTAC, Review 110, v. 11, 1957, p. 274. Some of the results have mean- 


while appeared in “Ona generalization of the prime pair problem,” by Herschel F. Smith, MTAC, 
v. 11, 1957, p. 274. 
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105 1091265 2839935 6503595 8741145 
16065 1615845 3243345 7187775 
19425 1954365 3400215 7641375 
43785 2822715 6005895 8062005 


Seven numbers out of 21 can be prime, eight out of 27 and nine out of 31. 
These all occur and are listed below. It is possible for ten numbers out of 33 to 
be prime and eleven out of 37, but these do not occur in the range; included in 
the list are such groups of nine out of 33 and ten out of up to 37 as occur in the 
range. 


[ane aw, Eh Lae ten me Ss St. SS 2D 
out of 
1271 31 19 — -_ F*— *——*F—-——- F*— * F— * JF — ZY 
5621 7 *— 17 13 —43——**—— *— * *— * * 21 
88781 747— 19 *— *——**—— *— * *— * * 24 27 31 
113141 7 * — -_ Fe— Fee FFT —- lO TF CU Fh FT * 11 — 27 1 FS 
165701 * * — _ *— *— — * * — — 103 — 53 17 — 11 , = 
284711 711 — 23101 — *——**—— *— * *— * * 24 27 
416381 7 11 — -_ F*— *—— FF —-— *— 2H F— * * — — 33 
626591 7 11 — *- *—17——**——- *— * *—~— * * 24 — 33 
855701 719— 149 *— *——**——- *— * *— * * 21 27 31 
1068701 * * — *— *——**—— 11 — 1947 — 29 7 21 
1146761 7 *— 43 13 — *——**—— *— * *— * * 24 27 
2580641 7 13 — -_ F*— *——*F*F—-—- *— * * — 11 457 — 27 
6560981 7 11 — * 290— *——**——- *— * *— * * 21 27 33 
6937931 7 * — ._ + *— FF 43—- * *—-— * * — — 33 37 
7540421 7 11 — 1879 13 —41——**—— *— * *— * * 21 
8573411 7 *— 191 *—17——**—— *— * *— * * 21 


In each line, the number in the left hand column is the first number for the 
line, an asterisk indicates a prime, a dash a multiple of 3 or 5 and a number the 
least factor of any other composite. There are thus seen to be eleven groups of 
seven primes in 21 numbers, eight of eight primes in 27, five of nine primes in 31 
(two in the line beginning with 113141) and four more of nine in 33, and one group 
of ten in 35 and one of ten in 37. The range 1 to 50 is excluded as being altogether 
exceptional. A list of the groups of four and five has been deposited in the VU MT 
file of MTAC (see Review 110, MTAC, v. 11, 1957, p. 274). 


King’s College 
Cambridge, England 
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60(A, B, C, D, H,I, K, L, M, P, R, V, X, Z_].—Georce E. Forsytue, Bibliography 
of Russian Mathematics Books, Chelsea Publishing Co., New York, 1956, 
106 p., 20 cm. Price $3.95. 


Professor Forsythe precedes his very useful bibliography of Russian mathe- 
matical books with an informative introduction which contains, among others, 
a complete review of the book. To quote Prof. Forsythe: 

“The subject matter of the books listed is mathematics, pure and applied, 
including tables beyond the most elementary, but excluding descriptive geometry. 
There are a few titles on quantum mechanics and other branches of mathematical 
physics, and more on mechanics, mathematical machines and nomography, but 
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these topics are far from completely covered. For textbooks, the subject matter 
is more advanced than calculus. 

The list is confined to books, excluding dissertations. Because of their mathe- 
matical interest, most of the serial monographs (Trudy) of the Steklov Mathe- 
matical Institute have also been included. The books in the list have been 
published (or reprinted) in Russian or Ukrainian since 1930. Translations into 
Russian are omitted.” 

The books listed cover a wide range of topics and are in very many cases 
written by some of the very best research men in Russia. Thus there is every 
reason for supporting Prof. Forsythe’s thesis that these books ‘‘would be of the 
greatest value to mathematicians and students.” If, as Prof. Forsythe seems to 
think, ignorance of the existence of these fine (and cheap) books and ignorance of 
how to obtain them are the major obstructions to their popularity, then the 
present bibliography ought to resuit in a run on Soviet mathematical books 
and periodicals. 

ABRAHAM SHENITZER 


Rutgers University 
New Brunswick, New Jersey 


61[A, B, C, D, E, F, K, L, M, N].—CRC Standard Mathematical Tables, Tenth & 
Eleventh Editions (Formerly Mathematical Tables from Handbook of Chem- 
istry and Physics), Chemical Rubber Publishing Company, Cleveland, Ohio, 
1957, ix + 480 p., 19.5 cm. Price $3.00. 


The section on Mathematical Tables from the Handbook of Chemistry and 
Physics considerably enlarged and greatly improved both in content and in the 
format of its tables made its first appearance in January 1954 as the Tenth Edi- 
tion of the C. R. C. Standard Mathematical Tables. 

Among the major improvements in the Tenth Edition are the following: 


1. A new section on vector analysis 
2. A table of Laplace transforms 
3. A table of Bessel Functions of orders zero and one, Hyperbolic Bessel 
Functions and Spherical Bessel Functions 
4. A table of sine, cosine and exponential integrals 
5. A table of Elliptic Integrals of the First and Second Kinds 
6. Extension of the tables of the hyperbolic functions sinh x, cosh x, and 
tanh x tox = 10 
7. Extension of the table of the exponential function to x = 10 
8. The revision and enlargement of the table of integrals to include an addi- 
tional ninety-four formulas 
9. A table of x? 
10. A table of F and ¢ for 1% and 5% distributions 
11. Commissioners 1941 Standard Ordinary Mortality Table with its auxiliary 
table of the commutation symbols at 23%. 
In addition to the above, the following are included in the Eleventh Edition: 
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1. An enlarged section on Fourier series with special power series, Bernoulli 
numbers and finite sine and cosine transforms 

. A new section on partial fractions with examples 

. A table of multiples of 7/2 

. A table of sums of the powers of natural numbers 

. A method for the tenfold extension of the tables of Factors and Primes 

. A table of particular solutions to ordinary and partial differential equations 
with constant coefficients. 


nan > wd 


These and other minor changes tend to make this volume more legible and 
more serviceable to all. 
I. A. STEGUN 


National Bureau of Standards 
Washington, D. C. 


62[B].—Toxyo NuMERICAL COMPUTATION BuREAU, Report No. 10. Table of 
Square Roots of Complex Numbers, 1956, i + 23 p., 25 cm. Not for sale. 


This table, with explanation by T. Sasaki, lists the real and imaginary parts 
of the square roots of 1 + ix and x + ito 11D for x = 0(.002)1. Second differences 
are given. Legibility is only moderate, but the table may be useful, because 
existing tables for the range 0 < x < 1 give fewer arguments and far fewer 
decimals; see Fletcher, Miller, and Rosenhead Index [1]. 

A. F. 
1. A. FLetcuer, J. C. P. Mrtier, & L. RosenneaD, An Index of Mathematical Tables, Scien- 


tte Copenting Service Limited, London, 1946, p. 33 (MTAC, Review 233, v. 2, 1946-47, 
p. 13-18). 


63[D, L].—G. E. Rreyno.ps, Table of (sin x)/x, Antenna Laboratory, Electronics 
Research Directorate, Air Force Cambridge Research Center, Bedford, Mass., 
AFCRC-TR-57-103, ASTIA Document No. AD117063, March 1957, iii + 5 
+ 200 p., 26 cm. 


This work was prompted by the frequent appearance of the function 
(sin x)/x in the mathematics of antenna design. Values are given to 8D for 
x = 0(.001)49.999. The table is machine printed in block form, the digit in the 
third decimal place of x being given as the top argument, and the earlier digits 
as left (and right) argument. There are no differences. The table was printed from 
punched cards produced by Datamatic Corporation of Newton Highlands, Mass. 
and the AFCRC Statistical Services Division. 

As far as the reviewer’s information goes, there are only two other tables 
which, although slighter on the whole, complement this new work to a significant 
extent. A manuscript table by C. Blanc [1] gives (sin x)/x, as well as its even 
derivatives to order 16, to 10D for x = 0(.01)4, and is worth mention in the 
present context because of its ten decimal places. A well-known volume by 
K. Hayashi [2] includes a table of (sin x)/x to 8D for x = 0(.01)10(.1)20(1)100, 
and is worth mention because of its values for x = 50(1)100 and its general rele- 
vance. Hayashi’s table, which contains errors, is very suitable for comparison 
with the new table, since both are strictly eight-place. The reviewer compared 
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some scores of corrected Hayashi values without finding any error in those of 
Reynolds. 


A. F. 
1. C. BLanc, see MTAC, v. 7, 1953, p. 188 (UMT 160). 
2. K. Hayvasui, Tafelu der Besselschen, Theta-, Kugel- und anderer Funktionen, Berlin, 1930, 
p. 30. See MTAC, v. 1, 1943-45, Review 12, p. 4. 


64[(D ].—CLovis FaucuEr, Tables trigonometriques contenant les valeurs naturelles 
des sinus et des cosinus de centigrade en centigrade du quadrant avec dix decimales, 
Gauthier-Villars, Paris, 1957, 51 p., 27 cm. Price 500 francs. 


These tables, of sin x and cosx to 10D for x = 0(0*.01)50® with first differ- 
ences, have been photo-offset from manuscript. The author states (in effect) that 
the values were derived by chopping the last two figures of 12D values calculated 
by quadratic interpolation in the fundamental tables of Andoyer. A dot appended 
after the 10th digit indicates that the value should be rounded up. The only 
checking specifically mentioned is a comparison with an existing 8D table. These 
tables, so far as is known to the reviewer, are the only published ones to 10D 
with argument interval 0*.01. (Peters has such a table in manuscript form.) Linear 
interpolation gives 8D, but full accuracy is attainable only with quadratic 
interpolation. 

T. H. SouTHARD 
University of California 
Los Angeles, California 


65(H, X].—JEAN PELTIER, Resolution Numerique des Equations Algebriques, 
Gauthier-Villars, Paris, 1957, iv + 244 p., 21 cm. Price $7.29. 


The equations considered are polynomial equations of high degree in a single 
variable, and the treatment is from the standpoint of the desk-machine user. This 
means that to some extent the book is behind the times because it is usually 
quicker and cheaper to mail such problems to an organization equipped with a 
high-speed automatic computer and set of subroutines for iterative root-finding. 

Only the most powerful methods are described and the reviewer is in agree- 
ment with the choice made. The opening chapters deal with elementary algebraic 
properties of the roots and the arithmetic operations required in the evaluation, 
multiplication and division of polynomials. Chapters 3 and 4 are devoted to the 
calculation of the moduli of the roots by Graeffe’s root-squaring process and the 
subsequent determination of the phases by the ‘‘highest common factor’’ method. 
An omission here is the use of the latter method merely to resclve the finite 
number of ambiguities in the determination of the accurate phases from the root- 
squaring computations. Chapter 5 deals with iterative methods, and later chapters 
cover error propagaton and equations with complex coefficients. The recommended 
general computational procedure is summarized in the final chapter. 

Many numerical examples are included, but little indication is given of check- 
ing procedures. In particular, the use of current checks on the root-squaring 
computations, which is vital in hand computing, does not appear to be mentioned. 


F. W. J. OLVER 


National Bureau of Standards 
Washington 25, D. C. 
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66[1].—L. N. Karmazina & L. V. Kurochxina, Tabliisy interpoltaisionnykh 
koéfflisientov (Tables of interpolation coefficients), Press of the Academy of 
Sciences of the USSR, Moscow, 1956, 376 p., 26.5 cm. Price 37.65 rubles. 


These tables have been prepared at the Computing Centre of the Academy 
of Sciences. The coefficients are given throughout to 10 decimals. The fractional 
part of the interval being , the following tables of Lagrangean coefficients are 
given: Table I, 3-point, p = —1(0.001)1; II, 4-point, 6 = —1(0.001)2; III, 
5-point, p = —2(0.001)2; IV, 6-point, p = —2(0.01)0(0.001)1(0.01)3; V, 7-point, 
p = —3(0.01) —1(0.001)1(0.01)3; VI, 8-point, = —3(0.01)0(0.001)1 (0.01)4. 
Table VII gives 3 to 11-point Lagrangean coefficients for sub-tabulation to tenths. 

All the coefficients for interpolation by differences are at interval 0.001, except 
Bessel’s which are at 0.0001 ; coefficients up to those of the seventh difference are 
given throughout and the error is said not to exceed 6 X 10-". The even Bessel 
coefficients are those of mean, not double mean, differences (i.e. they are 2B, 
2B,, 2B, in the notation of Interpolation and Allied Tables, London, 1956). 
Table VIII gives the Newton-Gregory coefficients, [X gives Bessel’s, X gives 
Stirling’s (—0.5 < p < 0.5), XI gives Everett’s (to Es). None of the tables of 
coefficients is provided with differences. 

All the coefficients have been newly computed except the Lagrangean ones 
which, apart from the portion of the eight-point table at interval 0.01, have been 
taken from the National Bureau of Standards Tables (1944). One may hazard a 
guess that the Russian equivalent of the National acco inting machines has eight 
registers; it is hard to see what else should make one stop at B;; the predilection 
for Bessel’s formula, which agrees with this reviewer’s personal preference, is 
noteworthy; more striking is the absence from the Introduction of any mention 
of throw-back. In view of Chebyshev’s connection with it, this is strange indeed. 


J. C. E. JENNINGS 
Birbeck College 
University of London 
London, W.C. 1., England 


67[1, X].—H. E. Sauzer, “Formulas for calculating Fourier coefficients,’’ Jn. 
Math. and Physics, v. 36, 1957, pp. 96-98. 


If a function is defined in the range (0, 27) as the polynomial of degree not 
exceeding eight which takes nine given values f, at the nine equally spaced argu- 
ments x = nx/4, where n = 0(1)8, then all the coefficients in the Fourier expan- 
sion of the function in the range (0, 27) are defined, and may be calculated from 
the given values by using equations of the form 


a1(fo + fs) + a2(fi + fr) + as(fe + fe) 
+ aa(fs + fs) + asf, 
bi(fo — fs) + b2(fi — fr) + b3(fe — fe) + ba(fs — fs). 


1 Qe 
— f f(x) cos nxdx 
2 Jo 


1 Qe 
— f f(x) sin nxdx 
2x Jo 


The table on page 97 lists values of the nine coefficients a; and 6; to 5D for 
n = 0(1)24, so that means are provided for calculating all Fourier coefficients up 
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to those: of cos 24x and sin 24x inclusive. The values of a; and b; tabulated are 
said to be “most probably correct to within around 1} units in the 5th decimal, 
even though they cannot be absolutely guaranteed to less than 3} units in the 
5th decimal.” 


A. F. 


68(I).—H. E. Satzer & PEGcy T. RoBerson, Table of Coefficients for Obtaining 
the Second Derivative without Differences, Convari-Astronautics, San Diego 12, 
California, 1957, ix + 25 p., 26 cm., deposited in thel/MT files. 


If m function values are given at m equidistant arguments at interval h, this 
table lists coefficients for determining the value of the second derivative at any 
point in the range of width (m — 1)h. The tables relate to the cases m = 5(1)9, for 
which the Lagrange interpolation polynomials are of degrees 4(1)8 respectively, 
and their second derivatives consequently of degrees 2(1)6 respectively. For 
n = 5(1)7, the coefficients are given for hundredths of h; for n = 8, 9, they are 
given for tenths of 4. The coefficients are in each case multiplied by an integer 
(6, 12, 360, 360, 10080 for m = 5(1)9 respectiveiy) such that the tabulated values, 
instead of involving recurring decimals, are exact in a convenient number of 
decimals (4, 6, 8, 5, 6 respectively.) This useful table was calculated on the IBM 
650 Magnetic Drum Data-Processing Machine. 

A. F. 


69[1].—H. O. Rosay, Interpolation Coefficients (=) for Newton’s Binomial 


Interpolation Formula, 20 tabulated sheets, 28 XK 38 cm., deposited in the 
UMT Fre. - 


Binomial coefficients (=) ,K = 1(1)5, S = 0(.01)5, 11D; however in many 


cases the last three digits are not reliable and accuracy is probably 8D. 

The tables were computed on SWAC for use with the Gregory-Newton inter- 
polation formula ([1], p. 96, equation 4.3.11). Similar tables for K = 1(1)6, 
S = 0(.01)100, 7D are in [2] according to [3]. 

T. H. SouTHarpD 


University of California 
Los Angeles, California 


1. F. B. HILDEBRAND, Introduction to Numerical Analysis, McGraw-Hill, New York, 1956. 
2. G. Veca & J. A. Hittsse, Sammlung Mathematischer Tafeln, Weidmann, Leipzig, 1933. 


3. A. V. LeBepev & R. M. Feporova, Spravochnik po Matematicheskim Tablitsam, Moscow, 
1956. [MTAC, Rev 49, v. 11, 1957, p. 104-106.] 


70[L].—A. A. Abramov, Tabliisy In T(z) v Kompleksnot oblasti (Tables of In T(z) 
in a complex region), Akad. Nauk SSSR, Moscow, 1953, 333 p., 26 cm. 
Price $3.50. 


This volume gives 6D tables of the real and imaginary parts of In T(x + zy) 
with second differences for x = 1(.01)2 and y = 0(.01)4. Real and imaginary 
parts are tabulated on opposite pages, each page having six columns corresponding 
to six values of x (the columns x = 1.05(.05)1.95 are repeated occurring as last 
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columns on one page and first columns on the following page) and 51 lines corre- 
sponding to fixed y (with repetitions for y = .5(.5)3.5). The introduction contains 
formulas for computation, error estimates, references both to books and tables 
(not including, however, Sibagaki’s tables mentioned below), and a nomogram 
for 3(£ — n)(& + 9 — 1)As. The volume is well printed. 

For other tables of the gamma function in the complex domain see RMT 234 
(MTAC, v. 2, p. 19), RMT 855 (MTAC, v. 5, p. 25-26), and RMT 1143 (MTAC, 
v. 7, p. 246). A small number of obvious misprints in the present work are cor- 
rected in pencil in the reviewer’s copy. Comparison of 1300 values given both in 
the present tables and those by Sibagaki (RMT 1143) showed about 100 dis- 
crepancies, mostly of one unit of the last decimal. 

A. Erpétyi 


California Institute of Technology 
Pasadena, California 


71[L ].—Car_-Ertk FrOBERG, Complete Elliptic Integrals, Lund University De- 
partment of Numerical Analysis, Table No. 2, G. W. K. Gleerup, Lund, 
Sweden, 1957, 82 p., 25 cm. Price 12 Kr. 


Ten decimal tables of the complete elliptic integrals of the first and second 
kinds as functions of the modulus. The functions tabulated are K(k), E(k), 
K'(k) = K(k’), E’(k) = E(k’) and two auxiliary functions 


S(k) = K'(k) — Ing, S’'(k) = K(k) — In = 
The intervals are k = 0.000(0.001) 0.900, 0.9000(0.0001) 1.0000. Second differ- 
ences, modified where necessary, are given for the ranges in which Everett's 
second order interpolation formula will give full accuracy. Elsewhere differences 
are omitted and auxiliary formulae are given, for use at the extreme ends of the 
table. In point of fact these auxiliary formulae could have been dispensed with 
had the author observed that while his function S(k) is not interpolable at the 
given interval for k < 0.005, the function 


K'(k) — =r In ~ 

k 
is interpolable without restriction. Since the function S(k) fails just where it is 
most needed, it is difficult to see the point of introducing it throughout the range 
of k, or indeed at all. 

The tables were calculated in two ways namely by Landen’s transformation 
and by power series but no information is given of the comparative labor of the 
two methods nor of the agreement of the results. 

The argument which presents itself naturally in calculations with elliptic 
functions and integrals is k® rather than k. Furthermore the use of k® enables the 
size of the table to be halved without loss. 

L. M. MILNE-THOMSON 


Brown University 
Providence, Rhode Island 
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72(L ].—Cari-Erik FrOBERG & Hans WILHELMssON, “Table of the function 
F(a, b) = Sf o* Ji(x) (x? + B*)-tdx,”” Fysiografiska Sallskapets I Lund Fér- 
handlingar, Bd. 27, 1957, p. 201-215. 

The function F(a, 5) defined in the title is tabulated to 6D (which for most of 
the given values means 6S) for a = .1(.1)2(.2)10 and b = 0(.1)2(.2)10. 

A short introduction mentions some physical applications of the function 
F(a, 6) and lists some limiting and special cases. Five-point Lagrangian inter- 
polation is said to give full accuracy except for Max (a, b) < 1; to facilitate 
interpolation ia this region, an auxiliary function f(a, b) defined by 


F(a, b) = 3L@ + 0) — b] — [f@, 5) 


is tabulated for a = .1(.1)1, 6 = 0(.1)1, 6D. 

The tables were calculated on SMIL, the electronic computer of Lund Uni- 
versity, by numerical quadrature. For 6 < 2 a power series expansion was used 
as a check. The error is said to be in general less than 1 unit in the last place. 


PETER HENRICI 
University of California 
Los Angeles, California 


73[L].—J. F. Nicnoras, A Table of fo7 exp (—1/u)du for Small Values of x. 
28 Mimeographed pages, 22 X 16.5 cm., deposited in the UMT File. 


Let I(x) = fo7 edu and f(x) = x~e!/*I(x). The table lists f(x) and I(x), 
x = .01(.0005).05075, 6S. The values given for f(x) are said to be accurate except 
for small errors in the sixth place, the values given for J(x) may be accurate to 
no more than four places for some values of the argument. 

The tables are intended for use in connection with chemical and other reac- 
tions which are thermally activated and governed by an equation of the form 
b = g(pje OF, 

They are related to other tables listed in the bibliography as well as a table 
whose existence is mentioned in a comment by Hastings [1 ]. They are also related 
to Hasting’s rational approximations for the function — Ei(—~) listed in [1 ]. This 
relation follows from the identity listed by the author, I(x) = xe~/* + Ei-/*, 

The tables are reproduced photographically from machine tabulations, with 
a 7 page introduction reproduced photographically from typescript. A few copies 
may be available directly from the author at The Division of Tribophysics, 
Commonwealth Scientific and Industrial Research Organization, University of 
Melborne, Australia. 

There is a minor misprint on page 3 of the introduction, where the author 
should have written I(x) = xe~/* — fo" ue" “du; he has a + sign between the 
two main terms. 

..' aoe 


1. Cecm HastinGs, Jr., JEANNE T. HAywarp, & JAMEs P. WonG, JR., Approximations for 
Digital Computers, Princeton, 1955, p. 187-190. See MTAC, v. 9, 1955, Review No. 56, p. 121-113. 
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74[L, P].—Rosert L. STERNBERG, J. S. Surpman, & S. R. Zonn, Table of the 
Bennett Functions A,w(h) and Amn(h,k), Laboratory for Electronics, Inc., 
Boston 14, Massachusetts, 18 p., 28 cm., deposited in U MT files. 


The announced purpose of these undated and unpublished tables is to provide 
for the simple or multiple Fourier series expansion of the output from a cut-off 
power law rectifier when responding to a simple or multiple frequency oscilla- 
tory input. 

The definitions of the functions are: 


2 
A,n(h) = “ x (cos u — h) cos mu du 


where 81: cosu>h, O0O<u<er 
2 
Amn(h, k) = SS. (cos u + k cosv — h) cos mu cos nv du dv, 


where R2:: cosu+kcosu>h, O0O<u, v<rz. 


Reference is made to two papers by W. R .Bennett [1]. 

A,,(h) is given to 6D for m = 0(1)9 and hk = 0(.05)1. A,na(h, R) is given to 
6D for m, m non-negative integers on the range 0 = m + m = 4andh = 0(0.1)2, 
k = 0.1(0.1)1. 

Computation was done on the IBM C.P.C. by from two to five different 
methods. All results are stated to have an absolute error of about 10~* at most in 
absolute value. 

Tuomas H. SouTHARD 


University of California 
Los Angeles, California 


1. W. R. Bennett, “New results in the calculation of modulation products,” Bell System 
Technical Journal, v. 12, 1933, pp. 228-243; also ‘The biased ideal rectifier," Bell System Technical 
Journal, v. 26, 1947, pp. 139-169. 


75(L, P].—Cu1ao-Min Cuu, Georce C. CLark, & Stuart W. CHURCHILL, 
Tables of Angular Distribution Coefficients for Light-Scattering by Spheres, 
University of Michigan Press, 1957, xv + 58 p., 22 X 28 cm. Price $3.00. 
Expressions for the differential scattering of light by a spherical particle are 

well known. Quite recently the present authors, following a suggestion of W. 

Hartel, derived forms which facilitate both numerical evaluation and application. 

The normalized differential scattering cross section is expressed by an expansion 

in Legendre polynomials. The coefficients are the quantities of interest here. 

They depend on two parameters, one the index of refraction of the medium 

surrounding the particle, the other a quantity, a, expressing the ratio of particle 

diameter to wave length. 

The coefficients are evaluated for a series of fifteen refractive indices, ranging 
from 0.90 to 2.00 and also for the limit at infinity. The values of a are selected 
integral ones from the range 1-30. The accuracy is usually +1 in the fifth decimal 
digit. The output of the computer is directly printed without transcription. It is 
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not explicitly stated to what extent the computations have been checked ; such 
checks should be straightforward here. Perhaps Part II should have been labelled 
“‘Total Scattering Coefficients’ rather than ‘Total Scattering Cross Sections.” It 
would also have been interesting to include some of the computing times. 

Together with the table of Legendre polynomials prepared by two of the au- 
thors (Clark and Churchill), calculations of differential scattering of light should 
now be simple and straightforward. 

N. METROPOLIS 

University of Chicago 
Chicago, Ill. 


76[M, P].—Gustav Doetscu, Anleitung zum praktischen Gebrauch der Laplace- 

Transformation, R. Oldenbourg Verlag, Munich 8, Germany, 1956, 198 p., 

23 cm. Price DM 22.00. 

The author of this work has written several well-known books on the Laplace 
transformation and its applications. The present book was undertaken at the 
request of the publishers and aims at providing a handbook of the Laplace trans- 
formation for the use of engineers, special consideration being given to the needs 
of the theory of servomechanisms. The presentation is lucid, and the style is 
free of mathematical sophistication ; at the same time the statements are precise 
and mathematically sound. Proofs are not given, except where they contribute 
to a full understanding of the results. Conditions of validity are not only stated 
but also fully discussed; and the reader is often warned against neglecting to 
observe proper caution in using these methods. 

After introducing the Laplace transformation in chapter 1, the so-called rules 
are discussed in chapter 2, applications to ordinary linear differential equations 
with constant coefficients, and to systems of such equations, are given in chapter 3, 
difference equations and recurrence relations are presented in chapter 4, partial 
differential equations in chapter 5, and integral equations and integral relations 
in chapter 6. The problem of inversion is discussed in chapter 7, and asymptotic 
expansions and stability in chapter 8. 

An Appendix of some sixty pages, almost one third of the book, was compiled 
by Rudolf Herschel and gives transform pairs with emphasis on those needed in 
servoengineering and other engineering problems. There are 41 general formulas 
followed by nearly 300 transform pairs arranged according to image functions 
(116 rational functions, 109 algebraic and elementary transcendental functions, 
and 56 functions arising in the solution of differential equations). 

Since the theoretical part covers a wide field in something like 130 pages, it 
is inevitable that not all topics are covered thoroughly. The presentation of ordi- 
nary differential equations with constant coefficients is very satisfactory as are, 
for a book of this kind, the unusually thorough discussions of the inversion 
problem and of asymptotic expansions. By comparison, the chapter on partial 
differential equations is somewhat meagre, and that on integral equations, skimpy. 
Differential and integral equations with time-lag are not discussed at all. 

In the field which it does cover thoroughly, the book will prove a reliable guide 








~~ > -. -— —_- 








ch 


It 


u- 


ld 





REVIEWS AND DESCRIPTIONS OF TABLES AND BOOKS 155 


to the engineer whose mathematical training enables him to understand fully 
simple and precise mathematical statements, and who needs precise and un- 
ambiguous instructions for using Laplace transforms. 

A. Erpéty1 
California Institute of Technology 
Pasadena, California 


77[S ].—MANFRED VON ARDENNE, Tabellen der Elektronenphysik Ionenphysik und 
Ubermikroscopie, 1. Band. Hauptgebiete, Deutscher Verlag der Wissen- 
schaften, Berlin, 1956, xvi + 614 p., 27 cm. Price DM 74. 


These two volumes contain a very complete collection of formulas, tables, 
references, and selections from the very extensive physical literature on electron 
physics, ion physics and the application of these subjects to microscopy. The 
material contained in these books relating to mathematics and mathematical 
methods is very meager. The latter topic is confined to a brief discussion (2 pages) 
of the Laplace equation and a reference to the Liebmann method for integrating 
this equation numerically. 

Volume II contains six pages devoted to mathematical formulas and refer- 
ences. Of these, two pages are devoted to references, one page to a four place 
tables of exponential functions, two pages to formulas pertaining to statistics, 
one page to a figure illustrating the logarithmic scale used in the book, and one 
page giving the first few terms of the MacLaurin expansion of various simple 
functions. 

A. H. T. 


78[V, X].—B. Erxin, Numerical Integration Methods or Supersonic Wings in 
Steady and Oscillatory Motion, UTIA Report No. 36, Institute of Aerophysics, 
University of Toronto, 1955, v + 37 p. + 22 figures + 4 tables, 28 cm. 
Available only on an exchange basis. 

The paper deals with the problem of the determination of a perturbation 
velocity potential on supersonic wings with such a geometry that analytic methods 
are not applicable and, therefore, a numerical integration is unavoidable. It is 
based on the linearized theory of supersonic flow, and the cases of a supersonic 
as well as of a subsonic leading edge are considered. 

The first part contains the results of the author’s two previous papers on the 
problem of steady flow. The essential novelty consists of an extension of the 
method for oscillating wings. Insofar as the reduced frequency is concerned, three 
different stages of complexity appear. For a very slow oscillation the computation 
formulae are nearly identical with those for steady motion. If the frequency is 
sufficiently small to permit the linearization of the equations with respect to the 
frequency, then the formulae generated are still relatively simple. Finally, if this 
simplification is not possible, then the method still gives a system of equations, but 
the process for their numerical solution seems to be rather complicated. 

Emphasis is placed on the presentation of the method in a form suitable for 
automatic digital computation. In some simple cases, such as those involving a 
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flat or a cambered delta wing, the numerical results are compared with the exact 
(linear) solutions and the agreement is found to be very good. 


PENTTI LAASONEN 
Kasarmikatu 2B11 


Helsinki, Finland 


79C[W, Z].—RicHarD G. CANNING, Electronic Data Processing for Business and 
Industry, John Wiley and Sons Inc., New York, 1956, xi + 332 p., 23 cm. 
Price $7.00 


The objective of this book ‘“‘is the outline of a program of study and planning 
leading to the preparation of a proposal to top management.”’ As the germ of his 
method, the author advocates (sic) systems engineering which he interprets for 
management personnel. The book is readable, many examples are given and the 
individual parts fit together quite well. Furthermore, the author recognized that 
hardware aspects could be expected to change (as they have) and accordingly 
qualified his material. This reviewer believes that the author has met his objec- 
tive; this book probably will continue to prove useful to management personnel, 
particularly if it passes through revisions designed to keep it abreast of both 
technology and experience. 

This reviewer does not wish to criticize this book on the basis of what it was 
not designed to provide. He does believe that it represents a reasonable set of 
case studies for clerical systems problems. The book is thus recommended for 
examination with the hope that such an observation process may stimulate or 
assist work on underlying numerical techniques or ingenious gadgetry, or both. 


W. H. MarLow 
The George Washington University 
Logistics Research Project 


s0[W, Z].—RicHarp G. CANNING, Installing Electronic Data Processing Systems, 
John Wiley and Sons Inc., New York, 1957, vii + 193 p., 22 cm. Price $6.00. 


Most managers who are considering Electronic Data Processing Systems, 
EDPS, should find this book of considerable value and interest. It should be of 
particular interest to those managers who have decided to include an EDPS in 
their organizations. The book is written in nontechnical language and assumes 
that the reader (the MANAGER) is relatively unacquainted with electronic 
computers. However, familiarity with such basic concepts as files, unit records, 
fields and some knowledge of programming is presumed. 

The author contends throughout the book . . . ‘‘that top management should 
not remain so aloof from the details of electronic data processing. . . . It is time 
for management to recognize that EDP is a management tool and it is the re- 
sponsibility of management to learn the advantages and limitations of this new 
tool. Perhaps the most effective way to learn about these new machines is to 
learn how to program a machine. Or to put it slightly differently, it is questionable 
whether a person really understands electronic data processing until he knows 
how to program.” This reviewer concurs with the above statements. 

This book consists essentially of a narrative description of the typical problems 
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encountered (and how they are solved) by an organization which is planning the 
installation of an EDPS. The organization is hypothetical (yet realistic); the 
problems encountered are not hypothetical, they are real problems. The discussion 
begins at a point in time where management has concluded a feasibility and 
application study and is beginning to plan for the installation of an EDPS. The 
discussion extends through a point in time (some three years later) where the 
installed system is in operation and the conversion period is almost complete. 
(The conversion means the changing over from the old operational system to the 
new electronic system.) 

The EDPS organization (be it a section or branch or department) is first fitted 
in perspective, i.e., relative to the entire organizational structure. Management 
then begins plans for acquiring and training the personnel needed to staff the new 
EDPS section. Included are detailed job descriptions, qualifications and basis 
for selection for each of the various types of personnel required to staff the EDPS 
section. Approximately one half of the book discusses the personnel problems, 
i.e., acquisition, training, duties and performance. The case history reveals that 
the programming effort required should not be underestimated, viz., one might 
expect programming efforts to expend one man hour per one single address in- 
struction coded and checked. The last half of the book discusses the problems 
encountered in the physical installation of the EDPS. This includes both hard- 
ware and personnel facilities and terminates with a thorough discussion of the 
problems encountered in the conversion period. 

In summary, this book should be of considerable value to anyone involved in 
the planning, installation, or operation of an Electronic Data Processing System. 


M. J. ROMANELLI 
Computing Laboratory 
Ballistic Research Laboratories 
Aberdeen Proving Ground, Maryland 


81[W, Z].—G. Kozmetsxky & P. Kircuer, Electronic Computers and Management 
Control, McGraw-Hill, New York, 1956, 296 p., 23 cm. Price $5.00. 


While it is difficult for this reviewer to judge, it would seem that this book, 
written primarily for the business executive and presupposing no technical train- 
ing, succeeds in its purpose ‘‘to explain certain new developments which may have 
a greater influence on the management of enterprise than any other single group 
of events have had since the first industrial revolution.” 

A principal shortcoming is a complete lack of documentation. There are four 
footnotes and at least a dozen quotations from various pronouncements and 
articles; however, this reviewer did not see a single precise bibliographical refer- 
ence within the text. In his opinion this greatly detracts from what seems to be 
a large collection of accurate and suggestive accounts of management experiences 
with computing equipments. So far as this reviewer could tell, the authors are 
reasonably meticulous in their limiting use of the present tense to computers and 
processes in being. 

The authors’ exposition on management control is capable and generally 
stimulating. Among the aspects described and related to computer applications 
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are the following: model-making, programming (planning), scheduling and feed- 
back, and simulation. 

In conclusion, this reviewer believes that the book would be of interest to 
those readers of this journal who wish to obtain some of the flavor of ‘‘business”’ 
data processing problems. He differs with the authors (p. 12) by believing that 
mathematicians must lead the way toward eventual resolution of difficulties in 
this particular area of numerical analysis. 

W. H. MarLtow 
The George Washington University 
Logistics Research Project 


82[W, X, Z ].—PurpvE UNIverRsITY COMPUTER RESEARCH PROGRAM, Proceedings 
of Symposium I, held at Purdue University, November 8 and 9, 1956, Purdue 
Research Foundation, iv + 76 p., 28 cm. 


The book is a transcript of seven papers presented at Purdue University 
Computer Research Program, a symposium held November 8 and 9, 1956; the 
table of contents follows: 

The Role of a University in an Industrial Society, by C. F. Kossack. 

Administration of Research, by R. A. Morgen. 

Reports on the Purdue Computer Research Program, by P. Brock. 

The Purdue Compiler, by Sylvia Orgel. 

Some Modern Linear Techniques in Practical Problems, by P. S. Dwyer. 

A Re-Evaluation of Computing Equipment Needs, by S. N. Alexander. 

Information Retrieval, by J. W. Mauchly. 

et 


83[X].—WiuiaM E. Mine, Numerical Solution of Differential Equations, John 
Wiley and Sons, Inc., New York, 1953, xi + 275 p., 23 cm. Price $7.25. 


The work under review has in the years since its appearance been widely used 
both as a text and as a reference book. The book owes this success to several very 
attractive features. Eminently teachable and lavishly illustrated with worked 
examples, it leads the student along a gentle path. While providing a great deal 
of useful information, it seldom exposes the reader to the stern discipline of ana- 
lytical rigor. Moreover, Professor Milne does not jump to the sterile conclusion 
that, because there is some good and some bad in every method, there is little to 
choose between them; on the contrary, he comes out with very definite opinions 
and thus offers firm guidance to the inexperienced worker in the field. 

If in the following we shall comment critically on some aspects of Milne’s 
book, we do not mean to belittle the author’s achievement, which was quite 
substantial at the time the book was written. Rather we believe that our remarks 
will illustrate the enormous development which has taken place in the field of 
numerical computation during the last decade and which Professor Milne was one 
of the first to help to bring about. 


Although high-speed machinery is mentioned at a few scattered places, it is 
probably fair to say that the whole outlook of Milne’s book is dominated by the 
idea that all computations are carried out on a desk machine. The faithful disciple 
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of the book is imagined to work doggedly with pencil and paper, erasing predicted 
values when replacing them by corrected ones. This basic philosophy of the book 
is, in this reviewer's opinion, a serious detriment to its value as a guide for use in 
connection with automatic computers. The following remarks may serve to 
substantiate this statement. 

The vexing problems that arise in the numerical integration of ordinary 
differential equations over a long range are not given adequate attention. The 
discussion of round-off is relegated to a rather superficial one-half page appendix. 
Although much attention is paid to the local error, no attempt is made to obtain 
realistic appraisals for the inherited error. Few numerical examples are carried 
further than ten steps. 

In the sections on ordinary differential equations, which together comprise 
about 100 pages, the emphasis is on difference methods. Methods of Runge-Kutta 
type are discussed on two pages (without proofs), and everything is done to 
discourage their use. The simple second order Runge-Kutta methods are not even 
mentioned. While it is true that difference methods are somewhat less involved 
from the point of view of pencil and paper work, this advantage is hardly relevant 
from the point of view of computing machinery. Here all genuine one-step methods 
(such as the various Runge-Kutta methods or their modern modifications) enjoy 
the tremendous advantage of being self-starting (and thus cutting the length of 
any code in half) and of rendering trivial the operation of changing the basic 
interval. 

‘Those who, in view of earlier publications of the author, expect the book to 
contain a wealth of relations involving finite differences will not be disappointed. 
Since for a given number of ordinates central difference formulas are the most 
accurate ones, these are particularly emphasized. Unfortunately these formulas 
are frequently unstable (in the sense of H. Rutishauser) if used for the integration 
of ordinary differential equations. A simple example may illustrate this. Milne 
discusses the classical Euler method (ya41 = yn + Hy’) in two lines. He then 
proceeds immediately to a discussion of the midpoint rule (ya41 = Ya-1 + 2hy’n) 
to which he devotes five pages. It so happens that the midpoint rule, quite apart 
from requiring a starting procedure, is unstable; it will not produce useful results 
for such a simple differential equation as y’ = —¥y over any reasonably long in- 
terval. Euler’s method, on the other hand, will produce good, if not very accurate 
results for any differential equation. The celebrated Milne method (known in 
European countries as method of central differences, or simply as Simpson’s rule) 
is also unstable, as was shown by Rutishauser [1 ]. While it is true that unstable 
methods usually produce good results over short ranges, their blind use in prob- 
lems involving many steps can lead to disastrous effects. 

The chapters on the solution of partial differential equations are handicapped 
by the fact that it was apparently not possible to give proper attention to several 
fundamental papers which became available shortly before the book was pub- 
lished. There is a long section on linear equations and matrices, intended to pro- 
vide the student with the tools for solving implicit partial difference equations in 
particular those arising from elliptic probiems. But the text does not mention 
the merits of the successive overrelaxation method, discussed independently by 
S. P. Frankel [2] and D. M. Young [3, 4], which has proved to be so eminently 
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successful in problems of this type. Again a great variety of partial difference 
operators is listed, but the discussion of numerical stability adds nothing to the 
1928 Courant-Friedrichs-Lewy paper. 

It has been pointed out that the general mathematical level of the book is 
elementary. To be specific, the level is well below the comparable book by L. 
Collatz [5]. While this is a desirable feature from the point of view of making 
the field accessible to students with a limited interest in mathematics, this 
reviewer has experienced a less desirable consequence of the lack of intellectually 
challenging material. To the bright mathematics student, if this is his first ac- 
quaintance with numerical analysis, the subject is made to appear a dull collection 
of recipes rather than a logical piece of mathematics. He has no opportunity to 
see the many interesting problems which merit his interest, and he will not feel 
attracted to a science which today more than ever needs his talents. 


PETER HENRICI 
University of California 
Los Angeles, California 


1. Von HEtnz RuTIsHAUSER, ‘Uber die instabilitat von mathoden zur integration gemohnlicher 
differentialgleichungen,” Zeit. f. ang. Math. u. Phys., v. 3, 1952, p. 65-74. 

2. STANLEY P. FRANKEL, “Convergence rates of iterative treatments of partial differential 
equations,” MTAC, v. 4, 1950, p. 65-75. 

3. D. M. Youn, Iterative Methods for Solving Partial Difference Equations of Elliptic Type, 
Doctoral thesis, Harvard University, 1950. 

4. D. M. Youns, “Iterative methods for solving partial difference equations of elliptic type,” 
Amer. Math. Soc., Trans., v. 76, 1954, p. 92-111. 

$. 4 re Numerische Behandlung von Differentialgleichungen, 2nd Ed., Springer-Verlag, 
Berlin, 1955. 


84[Z ].—Joun RoBERT Stock, Die mathematischen Grundlagen fiir die Organisation 
der elektronischen Rechenmaschine der Eidgenéssischen Technischen Hochschule, 
Inst. f. angew. Math., Miit., Zurich, No. 6, 1956, 73 p., 23 cm. Price Sw. 
Fr. 7.30. 


This seventy-three page book is the sixth in the current series on computers. 
and applied mathematics prepared by the Institute for Applied Mathematics at 
the Eidgenossishe Technische Hochschule in Zurich, Switzerland. 

The publication consists of a complete description of the ERMETH (Elek- 
tronische Rechenmaschine der Eidgenossishen Technischen Hochschule). In four 
Chapters, the mathematical and engineering characteristics, and the internal 
organization of the machine, are described. The first chapter presents a description 
of the mathematical characteristics, indicating number representation, instruction 
codes and machine operation. The second chapter gives the fundamentals of 
floating point, normalization and rounding, fixed point and index registers. The 
third chapter discusses the engineering fundamentals, including basic decision 
elements and arithmetic principles. The fourth chapter includes a discussion of 
internal organization and the manner in which the arithmetic unit performs 
addition, subtraction, multiplication and division. The control unit and error 
control methods are also described. 

A large volume of information is gi'ven in the five appendices. Appendix I 
contains a complete description of the external features: 
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TECHNICAL DATA: 


Drum storage capacity 
Word length 
Drum speed 


Pulse repetition rate in storage 

Pulse repetition rate in control and arithmetic unit 

Cycle time 

Average drum access time 

Average access time in the index—instruction 
counter—register 


OPERATION TIMES: 


Add, excluding storage access, average 
Multiply, excluding storage access, average 
Divide, excluding storage access, average 
Logical instruction 


COMPONENTS: 
Tubes 1,900 
Germanium diodes 7,000 


Of these, 40% are in the arithmetic unit 
25% are in the control unit 
20% are in the storage switching unit 
15% are in the input-output 


WORD STRUCTURE: 
Floating point number words 


Check digit (left end) 1 
Exponent 
Absolute value of mantissa (<10) 11 


Sign (right end) 


Fixed point number words 


Check digit (left end) 1 
Absolute value of mantissa (<1) 14 
Sign (right end) 


Instruction word 


Check digit (left end) 1 
Left or first instruction 7 
Operation 2 
Index 1 
Address 4 


10,000 words 
16 decimal digits 
6,000 RPM 
352 Ke 
32 Ke 
0.5 milliseconds 
5.0 milliseconds 


2.5 milliseconds 


4 milliseconds 
18 milliseconds 
28 milliseconds 

1 millisecond 
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Right or second instruction 7 
Operation 2 
Index 1 
Address 4 

Sigh (right end) 1 


Appendix II contains block diagrams, electrical schematics and logical dia- 
grams, of the arithmetic unit, including a description of the micro-steps performed 
by the arithmetic unit during the execution of various types of instructions. 

Appendix III gives similar information concerning the control unit. Included 
in Appendix IV is a detailed description of the networks contained in the adder. 

Appendix V contains four different typical programs that have been prepared 
for the computer. A bibliography is included, however all of the references are 
dated 1953 or earlier. 

The publication is written in clear, concise, easily understandable language. 
Although the ERMETH is a relatively slow computer (no high-speed storage 
other than the three arithmetic registers), the discussion and analysis of general 
computer organization and principles of operation are very interestingly pre- 
sented. The reader will be quite pleased with the manner in which such subjects 
as number systems, arithmetic schemes, instruction codes, and programming 
techniques are handled. Sufficient information is given to permit the reader to 
code problems for the ERMETH. The book is worthwhiie reading for anyone 
who wishes to familiarize himself with the general techniques and principles of 
computer construction, organization, and operation, even though the ERMETH 
itself may not be the reader’s particular concern. 

- Martin H. WEIK 
Computing Laboratories 


Ballistic Research Laboratories 
Aberdeen Proving Ground, Maryland 


85[Z ].—Taxkasui Kojima, The Japanese Abacus—Its Use and Theory, Charles 
E. Tuttle Co., Tokyo, Japan, 1957, 102 p., 18 cm. Price $1.25. 


There has been some interest in the abacus as the so-to-speak most elementary 
digital calculating instrument in use, but there has been no book in English which 
describes the abacus adequately. 

This book appears to be the first of that kind. 

Beginning with the story of the contest held in Japan in 1946 between oper- 
ators of the Japanese abacus and the electro-mechanical desk calculator, it deals 
briefly with the history of the abacus in general, then goes on to explain step- 
by-step manipulations of beads of the Japanese abacus to do four arithmetic 
operations. 

The book shows that by following various rules, the arithmetic operations 
(especially addition and subtraction) can be reduced to quite mechanical opera- 
tions of the abacus requiring very small amounts of mental effort (it is only 
necessary to remember 10’s complements of 1 — 9 and 5’s complements of 1 — 4). 

The book also contains some exercises for those who wish to acquire skill in 
manipulating the Japanese abacus. 
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The abacus is used quite extensively in Japanese business establishments. In 
the hands of experts, it is faster than the usual desk calculator for addition or 
subtraction and about as fast for multiplication or division of 10 to 12 digit 
numbers. The fact that the price of the Japanese abacus is approximately $3.00 
may be one of the big reasons for its popularity among the small business estab- 
lishments in Japan. However it must be borne in mind that to attain reasonable 
proficiency in the use of the abacus requires a number of months of practice. 

The monograph provides interesting reading for those who are curious about 
the Japanese abacus. 


MASANAO AOKI 
University of California 
Los Angeles, California 


TABLE ERRATA 


261.—Lre CENTRE NaTionaL D’Etupes pes TéLécommunications, Tables des 
fonctions de Legendre associées, Paris, 1952. [MTAC, v. 7, RMT 1110, p. 178; 
MTAC, v. 8, Table Erratum 233, p. 28.] 


The authors report the following errors: 


38, P%s,6(cos 17°) = —0.144118 and not —0.144072 
42, Ps 6(cos 17°) —0.261274 and not —0.261234 
82, P's, 6(cos 17°) —4,.053574 and not —4.053695 
86, P's, 6(cos 17°) = —3.47195 and not —3.47219 

p. 136, P%,6(cos 17°) = 49.29970 and not —49,30035. 


PUD? 


They further report that there are a number of instances of poor printing, 
which might lead to confusion. The list follows: 


p. 17, P;°(cos 80°) = —0.2473819—the 8 is illegible 

p. 18, P®; 9(cos 14°) = 0.7338195—the 9 is poorly printed and could be mis- 
taken tor an 8 

p. 60, P's 1(cos 13°) = —1.3328662—the 8 is poorly printed 

p. 78, P*z,9(cos 31°) = 2.068962—the 6 is poorly printed and could be mis- 
taken for an 8. 


Louis ROBIN 
Recherches Mathématiques 
L’Ingénieur en Chef des Télécommunications 
3 Avenue de la Republique 
Paris, France 


262.—CarL-ERIK FROBERG, Hexadecimal Conversion Tables, C. W. K. Gleerup, 
Lund, Sweden, 1957. [MTAC, Review 82, v. 11, 1957, p. 208; MTAC, Table 
Erratum, v. 11, 1957, p. 309. ] 


On p. 10, for 0.38 30A3D 70A3D 40A3D read 0.38 30A3D 70A3D 70A3D. 


B. ASKER 
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263.—D. R. KApREKAR, Cycles of Recurring Decimals, v. 1. (From N = 3 to 161 
and some other numbers.) Khare Wada, Deolali, India, 1950. Published by 
the author. [MTAC, Review 1126, v. 7, 1953, p. 238. ] 


A rather comprehensive examination of this table from N = 3-163, inclusive 
revealed a number of errors not listed in the author’s errata sheet at the back of 
the volume. Details of the errors noted by this investigator are given below. 

In addition, there seem to be errata in Kaprekar’s ‘‘errata’’ list, The correction 
he gives for page 16, line 20 should read 38-59-55 instead of 38-59-15; also the 
correction given for page 46, line 15 should read 115—93 instead of 115-98. The 
corrections given for page 8, line 8; also page 20, line 22 seem unnecessary ; these 
“‘corrections” had already been made in the volume seen. Kaprekar’s ‘‘correction”’ 
for page 24, line 13 is garbled; it should read 1, 1, 1, 8, 8, 8. 





Page Pair of Upper Line of Pair Lower Line of Pair 
of Lines in 
Volume N Cycle Cycle Reads Should Read Reads Should Read 
6 3S D9 iat 23-2 23-32 
37 5 ist 8, 6, 2 1, 6, 2 
15 we 2 1st Entire cycle wrong—See Note (a) 
16 107 2 2nd 8,9,0 9,9,0 
17 111 12 Ist 18-69-24 64-85-73 1,6,2 5,7, 6 (b) 
= wy FS te Entire cycle wrong—See Note (c) 
2.6123—CGtsié«éiAitst, 4,8,9 0,8,9 
22.6 61330—C 4st 48-15-17 68-15-17 
48 159 6 Ist 141-119-154 151-79-154 - 
8 ist 115-77-134 119-77-134 
49 163 2 ist 64-114-162 44-114-162 
2nd 122-179-138 122-79-138 
Notes 


(a) Cycle 25 improperly duplicates Cycle 22; it should read 45—46—56-55 


(b) Entire cycle wrong 
(c) Cycle 3 should read 31—76—-58-112-67-85 
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CHARLES R. SEXTON 


3009 Claremont Ave. 
Berkeley, Calif. 


264.—D. R. Kaprekar, Cycles of Recurring Decimals, v. Il. (From N = 167 to 
213 and many other numbers.) Khare Wada, Deolali, India, 1953. Published 
by the author. [MTAC, Review 1205, v. 8, 1954, p. 148.] 


In contrast to v. I, the v. II examined did not have an errata list. A rather 
comprehensive study of the tables from N = 167-213, inclusive, developed the 
errors listed below. 
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Page Pair of Upper Line of Pair Lower Line of Pair 
of Lines in 


Volume N Cycle Cycle Reads Should Read Reads Should Read 
1 167 1# 3rd =: 180-13-130 18-13~-130 





5th 4; %.1 ce | 
84-5-Blank 84-5-—50 5,0, blank 5,0,2 
S ik. 7 2nd 5, 8,8 9, 8,8 
& 7 % 2nd 3, 9,3 2,9, 3 
6 181 14## 5th 7, 6,4 4, 6,4 
: so 2nd 88-189-171 38-189-171 
193 1# 1st 9,5,1 0, 5,1 
3rd = 103-125-92 109-125-92 
11 201 3 3rd 5, 8,0 8, 8,0 
12 203 1 3rd > > 4, 3,3 
2 7th 143-39-187 146-39-187 
207 1 2nd 131-118-145 136-118-145 
| a: ae 2nd 9, 6, 3 7, 6, 3 
Notes 


# First half 
## Second half 


3009 Claremont Ave. 
Berkeley, Calif. 


CHARLES R. SEXTON 


NOTES 


National Bureau of Standards and National Science Foundation Training 
Program in Numerical Analysis for Senior University Staff 
Spring 1959 


The National Bureau of Standards is planning, conditional upon support by 
the National Science Foundation, to hold its second Training Program in Nu- 
merical Analysis for Senior University Staff. 
The purpose of this program is to give regular university staff a training in 
numerical analysis which will enable them to direct the operation of a university 
computing center, and to organize training and research in numerical analysis 
on their return to their own institutions. It will occupy the whole of the second 
semester of the academic year 1958-1959 (from February 9 to June 5) and has 
been arranged for that time so that participants may become familiar with the 
details of their own computing equipment during the following summer and be 
able to conduct courses in the academic year 1959-1960. 
Applications must be received not later than October 15, 1958. They should 
be addressed to: 
Dr. Philip J. Davis, Chief 
Numerical Analysis Section 
Applied Mathematics Division 
National Bureau of Standards 
Washington 25, D. C. 

and should include: 


1. The academic history of applicant and a list of his publications. 
2. A statement about the computational program of the institution. 





NOTES 


CORRIGENDA 


A. J. M. Hitcucock, ‘‘Polynomial approximations to bessel functions of 
order zero and one and to related functions,” MTAC, v. 11, 1957, p. 86, 
line 2 from bottom, 


for Jo(x) =1-— JE (Pace cos (« + *) — Qo(x) sin (« + *)) 


read Jo(x) =1— J2 (Po cos (« + =) — Qo(x) sin (« + *)) 








